Hello !
I am using for the first time the quadrature element. But I met an error with the NonlinearVariationalProblem I define due to the ‘Abs’ operator. I have already met this error when I used a linear solver but in this case I use a Newton solver. I define manually the Jacobian because it seems the derivature function does not work with Quadrature element.
I do a single egde notched test in mode II solved with a phase field model (staggered resolution).
I am facing the following error :
Form has no parts with arity 2.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File "/scratchm/lmersel/FEniCS/PFmodel/PF_fenics_Lamia/SingleNotchTest/Split_Amor/PetcsSNES/Hmod_QuadSpace/Main_SENT_SNES_Am_Hmod_QuadSpace.py", line 206, in <module>
problem_u = NonlinearVariationalProblem(f_u,u_,bcu , J_u) #==========> Problme
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfin/fem/problem.py", line 92, in __init__
J = Form(J, form_compiler_parameters=form_compiler_parameters)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfin/fem/form.py", line 43, in __init__
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
return local_jit(*args, **kwargs)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/jitcompiler.py", line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/jitcompiler.py", line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/compiler.py", line 142, in compile_form
return compile_ufl_objects(forms, "form", object_names,
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/compiler.py", line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
form_datas = tuple(_analyze_form(form, parameters)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/analysis.py", line 89, in <genexpr>
form_datas = tuple(_analyze_form(form, parameters)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ffc/analysis.py", line 169, in _analyze_form
form_data = compute_form_data(form,
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 84, in map_expr_dags
r = handlers[v._ufl_typecode_](v)
File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-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.
My code is :
#https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html#vonmisesplasticity
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
import sys
import time
level = 16
set_log_level(level)
ppos = lambda x: (x+abs(x))/2.
pneg = lambda x: (x-abs(x))/2.
def w(alpha):
return alpha**2
def w_prime(alpha):
return 2*alpha
#Degradation function
def gk(alpha, k= 1e-6):
return (1-alpha)**2 + k
def gk_prime(alpha, k= 1e-6):
return 2*(alpha-1)
def epsilon(u):
e = sym(nabla_grad(u))
return e
def psi_plus(u, K0,G):
eps = epsilon(u)
return 0.5* K0*ppos(tr(eps))**2 + G *inner(dev(eps),dev(eps))
def psi_moins(u, K0,G):
eps = epsilon(u)
return 0.5* K0*pneg(tr(eps))**2
def psi_0(u, G, K0):
return psi_plus(u, K0,G)+psi_moins(u, K0,G)
def psi(u,alpha, G, K0):
return gk(alpha, k= 1e-6)*psi_plus(u, K0,G)+psi_moins(u, K0,G)
def sigma_amor(u,alpha , G, K0,ndim):
eps = epsilon(u)
plus_tr_eps = ppos(tr(eps)) #0.5*(tr(eps)+abs(tr(eps)))
moins_tr_eps = pneg(tr(eps)) #0.5*(tr(eps)-abs(tr(eps)))
psi_plus_prime = 0.5* K0 * plus_tr_eps * Identity(ndim) + 2*G *dev(eps)
psi_moins_prime = 0.5* K0 * moins_tr_eps * Identity(ndim)
s = gk(alpha, k= 1e-6) * psi_plus_prime + psi_moins_prime
return s
# ------------------
# Parameters
# ------------------
ERROR = 40
set_log_level(ERROR) # log level
parameters.parse() # read paramaters from command line
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs" # quadrature is deprecated
"""
=========================================================
Physical parameters
=========================================================
"""
E, nu = Constant(210000), Constant(0.3) #[MPa]
G = Constant(E / (2.0*(1.0 + nu))) # [MPa]
lmbda = Constant (E * nu / ((1.0 + nu)* (1.0 - 2.0*nu)))
K0, Gc = Constant(lmbda+2./3.*G), Constant(2.7)
cw, lw = Constant(2), Constant(0.01) # [mm]
T = 0.02 # final time
num_steps = 200 # number of time steps
dt = T / num_steps # time step size
time_step = np.linspace(0, T, num_steps)
"""
=========================================================
Create mesh
=========================================================
"""
nx = 150
ny = 150 # 14 400 cell h = 0.01 mm
mesh = UnitSquareMesh(nx, ny,'crossed')
Nnode = len(mesh.coordinates())
Ncell = len(mesh.cells())
coord = mesh.coordinates()
ndim = mesh.topology().dim()
"""
=========================================================
Space discretization
=========================================================
"""
degU ,degDam = 1, 1
deg_history, deg_stress = 2, 2
V = VectorFunctionSpace(mesh,'Lagrange', degU)
V_alpha = FunctionSpace(mesh,'Lagrange', degDam)
#V_H = FunctionSpace(mesh,'DG', degH) #degH = 0
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_history, quad_scheme='default')
V_H = FunctionSpace(mesh, W0e)
We = TensorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, shape=(3,3), quad_scheme='default')
Tens = FunctionSpace(mesh, We)
"""
=========================================================
Boundary condition
=========================================================
"""
Disp_shear = Expression( ('t', '0.'), t = 0.0, degree = 1 ) # Constant((0.,0.))
# boolean boundary function
boundary_bottom = 'near(x[1], 0)'
boundary_up = 'near(x[1], 1)'
boundary_left = 'near(x[0], 0)'
boundary_right = 'near(x[0], 1)'
class Precrack(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (0., 0.5)) and near(x[1], 0.5) ) # and : &&
crack = Precrack()
bottom = CompiledSubDomain(boundary_bottom)
up = CompiledSubDomain(boundary_up)
left = CompiledSubDomain(boundary_left)
right = CompiledSubDomain(boundary_right)
boundaries = MeshFunction("size_t", mesh,1)
boundaries.set_all(0)
bottom.mark(boundaries, 1) # mark left as 1
up.mark(boundaries, 2) # mark right as 2
left.mark(boundaries, 3) # mark right as 3
right.mark(boundaries, 4) # mark right as 4
crack.mark(boundaries, 5) # mark right as 5
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(4), ds(0) all the edges
bcu_bottom = DirichletBC( V , Constant((0.,0.)), boundary_bottom)
bcu_right = DirichletBC( V.sub(1) , Constant(0.), boundary_right)
bcu_left = DirichletBC( V.sub(1) , Constant(0.), boundary_left)
bcu_up = DirichletBC( V , Disp_shear, boundary_up )
bcu = [bcu_bottom, bcu_up, bcu_right, bcu_left]
bcp_bottom = DirichletBC(V_alpha, Constant(0.), boundary_bottom)
bcp_up = DirichletBC(V_alpha, Constant(0.), boundary_up)
bcp_crack = DirichletBC(V_alpha, Constant(1.), boundaries,5)
bc_alpha = [bcp_bottom, bcp_up, bcp_crack]
"""
=========================================================
Define Trial and Test Function
Define functions for solutions at current times steps
=========================================================
"""
du,v,u_ = TrialFunction(V), TestFunction(V), Function(V)
dalpha,beta,alpha_ = TrialFunction(V_alpha), TestFunction(V_alpha), Function(V_alpha)
# INITIAL CONDITION
u_old, alpha_old = Function(V), Function(V_alpha)
H_new ,H_old = Function(V_H),Function(V_H)
g = gk(alpha_old,1.0e-6)
sig_funct = Function(Tens, name="Stress")
Psi_elastic = Function(V_H, name ="elastic_energy_density" )
"""
=========================================================
Define Variational problem
=========================================================
use Quadrature elements a custom integration measure must be defined to
match the quadrature degree and scheme used by the Quadrature elements
"""
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
#---------------- Mechanical equation
tn = Constant((0., 0.)) # external stress
f_u = inner( sigma_amor(u_, alpha_, G , K0, ndim) , epsilon(v) ) * dxm + dot( tn, v ) * ds
a_u = lhs(f_u)
L_u = rhs(f_u)
#J_u = derivative(f_u, u_, du) # jacobian
J_u = inner(epsilon(du), sigma_amor(v, alpha_, G , K0, ndim))*dxm # jacobian
problem_u = NonlinearVariationalProblem(f_u,u_,bcu , J_u)
solver_u = NonlinearVariationalSolver(problem_u)
newton_prm = solver_u.parameters['newton_solver']
newton_prm['relative_tolerance'] = 1e-8
newton_prm['absolute_tolerance'] = 1e-8
newton_prm['maximum_iterations'] = 2000
newton_prm['error_on_nonconvergence'] = False
#--------------- Phase Field equation
f_alpha = H_new*inner(beta, gk_prime(alpha_)) * dx \
+ (Gc/(cw*lw)) * (beta*w_prime(alpha_) + (2.*lw**2)*inner(grad(beta), grad(alpha_))) * dx
a_alpha = lhs(f_alpha)
L_alpha = rhs(f_alpha)
J_alpha = derivative(f_alpha, alpha_, dalpha) # jacobian
alpha_min = interpolate(Constant(DOLFIN_EPS), V_alpha) # lower bound
alpha_max = interpolate(Constant(1.0), V_alpha) # upper bound
problem_alpha = NonlinearVariationalProblem(f_alpha,alpha_ ,bc_alpha, J_alpha)
problem_alpha.set_bounds(alpha_min, alpha_max)
solver_alpha = NonlinearVariationalSolver(problem_alpha)
snes_prm = {"nonlinear_solver": "snes",
"snes_solver" : { "method": "vinewtonssls",
"maximum_iterations": 2000,
"relative_tolerance": 1e-9,
"absolute_tolerance": 1e-9,
"report": True,
"error_on_nonconvergence": False,
}}
solver_alpha.parameters.update(snes_prm)
# Stocked Forces and Energy value
forces = np.zeros((len(time_step),3))
loading = np.zeros((len(time_step),2))
energies = np.zeros((len(time_step),3))
#Data_SIG = np.zeros((num_steps, 3, Nnode))
Data_SIG = np.zeros((num_steps, 3, Ncell))
Data_Usol = np.zeros((num_steps, 2, Nnode))
H_tab = np.zeros((num_steps, Ncell))
#---------------- RESOLUTION
t= 0.
it = 0
tol = 0.001 # 1e-6
maxiter = 300
stop_tt = 10
import ufl
for t in time_step[:stop_tt]:
it +=1
Disp_shear.t = t
loading[it]= np.array([it,Disp_shear.t])
print("\n========================================================")
print("\n -> Timestep : %d s, Time : %.3g s, Ux_given : %.3g mm "%(it, t, Disp_shear.t))
print("\n========================================================")
k = 0
res_H = 1.
#A, Res = assemble_system(J_u, f_u, bcu)
#Du.interpolate(Constant((0, 0)))
print("\n ----- > Start Staggered Loop")
while (res_H > tol) and (k < maxiter):
k +=1
#----Step 1 : Displacement resolution
iter_u = solver_u.solve()
#solve(A, du.vector(), Res, "mumps")
#Du.assign(Du+du)
#----Step 2 : Computation of H energy history
elastic_energy_density=project(psi_plus(u_, K0,G), V_H)
H_new.vector().set_local(np.max((elastic_energy_density.vector().get_local(),H_tab[it-1]),axis=0))
print("\n________________________________________")
print("\n PETscSNESSolver Parameters for Damage pb")
print("__________________________________________")
iter_alpha = solver_alpha.solve()
print("\nEnd SNES Parameters__________________________")
#---- Check Residual
H_gap = H_new.vector().get_local() - H_old.vector().get_local()
index_max = np.argmax(H_gap)
H_gap_max = (H_new.vector().get_local() - H_old.vector().get_local()).max()
res_H = abs(H_gap_max)/ H_old.vector().get_local()[index_max]
err_u = errornorm(u_, u_old, norm_type = 'l2', mesh = None)
err_alpha = errornorm(alpha_, alpha_old, norm_type = 'l2', mesh = None)
err = max(err_u, err_alpha)
#---- Step 4 : Update previous solution
#u_old.assign(u_+Du)
u_old.assign(u_)
alpha_old.assign(alpha_)
g = gk(alpha_old,1.0e-6)
H_old.assign(H_new)
print( "\nIteration: %d, err_u: %.8g, err_alpha : %.8g, Residual: %.8g \n" % (k, err_u, err_alpha, res_H))
iter = [k]
print("\n ----- > Timestep : %d"%(it))
ux , uy = u_.sub(0, deepcopy=True), u_.sub(1, deepcopy=True)
H_tab[it, :] = project(H_new, V_H).vector().get_local()
Psi_elastic.assign(project(psi(u_,alpha_, G, K0), V_H))
# vertice value (mapping different from ddf value)
Data_Usol[it, 0, :] = ux.compute_vertex_values()[:]
Data_Usol[it, 1, :] = uy.compute_vertex_values()[:]
sig_funct.assign(project(sigma_amor(u_, alpha_, G , K0, ndim), Tens))
s11 = project(sigma_amor(u_, alpha_, G , K0, ndim)[0, 0], V_H)
Data_SIG[it, 0, :] = s11.vector().get_local()
s22 = project(sigma_amor(u_, alpha_, G , K0, ndim)[1, 1], V_H)
Data_SIG[it, 1, :] = s22.vector().get_local()
s12 = project(sigma_amor(u_, alpha_, G , K0, ndim)[1, 0], V_H)
Data_SIG[it, 2, :] = s12.vector().get_local()
# Calculate the axial and tangent force resultant
Fn = assemble(sigma_amor(u_,alpha_, G , K0, ndim)[0,0]*ds(2))
Ft = assemble(sigma_amor(u_,alpha_, G , K0, ndim)[1,1]*ds(2))
forces[it] = np.array([t,Fn, Ft])
#Calculate individually the energy
elastic_energy = 0.5*inner(sigma_amor(u_,alpha_, G , K0, ndim), epsilon(u_))*dxm
dissipated_energy = Gc/cw*(w(alpha_)/lw + lw*dot(grad(alpha_), grad(alpha_)))*dx
elastic_energy_value = assemble(elastic_energy)
surface_energy_value = assemble(dissipated_energy)
if it == 1 :
surf_off_set = surface_energy_value
print(' surf_offset_precrack : ', surf_off_set)
energies[it] = np.array([t,elastic_energy_value,surface_energy_value - surf_off_set])
# monitor the results
if MPI.comm_world.rank == 0:
print( "Eel_max: %.8g \nH_n_max: %.8g \nalpha_max: %.8g" % (elastic_energy_density.vector().max(), H_new.vector().max(), alpha_.vector().max()) )
print("\nIterationU and IterationAlpha: (%i,%i)"%(iter_u[0],iter_alpha[0]))
print("\nElastic and surface energies : (%g,%g)"%(elastic_energy_value,surface_energy_value - surf_off_set)) # WHY they display a scalar (a value) and not scale field (vector)
print("-------------------------------------------------------------")
Do you have an idea of how could I define the formulation differently to avoid the error ?
Thanks for the hints !
Lamia
