Hi all, recently I’m trying to learn some code to solve the thermoplastic problem with the following code:
import numpy as np
from arguments import args
from utils import walltime
import fenics as fe
import os
import glob
fe.parameters["form_compiler"]["representation"] = 'uflacs'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
@walltime
def simulation():
ambient_T = args.T_ambient
rho = args.rho
Cp = args.c_p
k = args.kappa_T
h = args.h_conv
eta = args.power_fraction
r = args.r_beam
P = args.power
EPS = 1e-8
E = fe.Constant(70e3)
nu = fe.Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = fe.Constant(250.)
Et = E/100. ### 切向模量
H = E*Et/(E-Et) ### 塑性模量，又称硬化系数
x0 = 0.2*args.domain_x
y0 = 0.5*args.domain_y
total_t = 1200*1e-6
vel = 0.6*args.domain_x/total_t
dt = 1e-6
ts = np.arange(0., total_t + dt, dt)
print(f"total time steps = {len(ts)}")
ele_size = 0.01
mesh = fe.BoxMesh(fe.Point(0., 0., 0.), fe.Point(args.domain_x, args.domain_y, args.domain_z),
round(args.domain_x/ele_size), round(args.domain_y/ele_size), round(args.domain_z/ele_size))
mesh_file = fe.File(f'data/vtk/pbf/mesh/mesh.pvd')
mesh_file << mesh
# Define bottom surface
class Bottom(fe.SubDomain):
def inside(self, x, on_boundary):
# The condition for a point x to be on bottom side is that x[2] < EPS
return on_boundary and x[2] < EPS
# Define top surface
class Top(fe.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[2] > args.domain_z - EPS
# Define the other four surfaces
class SurroundingSurfaces(fe.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[0] < EPS or x[0] > args.domain_x - EPS or x[1] < EPS or x[1] > args.domain_y - EPS)
bottom = Bottom()
top = Top()
surrounding_surfaces = SurroundingSurfaces()
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
bottom.mark(boundaries, 1)
top.mark(boundaries, 2)
surrounding_surfaces.mark(boundaries, 3)
ds = fe.Measure('ds')(subdomain_data=boundaries)
normal = fe.FacetNormal(mesh)
deg_u = 2
deg_stress = 2
# TODO: can we do a tensor element?
U = fe.VectorFunctionSpace(mesh, "CG", deg_u)
We = fe.VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=9, quad_scheme='default')
W = fe.FunctionSpace(mesh, We)
# We = fe.TensorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
# W = fe.FunctionSpace(mesh, We)
W0e = fe.FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = fe.FunctionSpace(mesh, W0e)
sig = fe.Function(W)
sig_old = fe.Function(W)
n_elas = fe.Function(W)
beta = fe.Function(W0)
# name = "a b" is bad (having a space)!
p = fe.Function(W0, name="Cumulative plastic strain")
u = fe.Function(U, name="Total displacement")
du = fe.Function(U, name="Iteration correction")
Du = fe.Function(U, name="Current increment")
v = fe.TrialFunction(U)
u_ = fe.TestFunction(U)
u_bcs = [fe.DirichletBC(U, fe.Constant((0., 0., 0.)), bottom)]
def eps(v):
e = fe.sym(fe.grad(v))
return e
def sigma(eps_el):
return lmbda*fe.tr(eps_el)*fe.Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return fe.as_tensor([[X[0], X[1], X[2]],
[X[3], X[4], X[5]],
[X[6], X[7], X[8]]])
def as_long_vector(X):
return fe.as_vector([X[0, 0], X[0, 1], X[0, 2], X[1, 0], X[1, 1], X[1, 2], X[2, 0], X[2, 1], X[2, 2]])
ppos = lambda x: (x+abs(x))/2.
def thermal_strain(dT):
alpha_V = 1e-5
return alpha_V*dT*fe.Identity(3)
def proj_sig(deps, dT, sig_old, p_old):
sig_n = as_3D_tensor(sig_old) ### 将旧的应力场转为三维张量
# sig_elas = sig_n + sigma(deps)
d_eps_T = thermal_strain(dT) ### 计算温度变化dT引起的热应变
sig_elas = sig_n + sigma(deps - d_eps_T) ### 计算当前弹性应力张量
s = fe.dev(sig_elas) ### 计算弹性应力张量的偏量部分
sig_eq = fe.sqrt(3/2.*fe.inner(s, s)) ### 计算等效应力
f_elas = sig_eq - sig0 - H*p_old ### 计算弹性力，用于评估是否需要激活塑性本构。>0,应力状态超过屈服应力；
dp = ppos(f_elas)/(3*mu+H) ### 计算塑性应变增量
n_elas = s/sig_eq*ppos(f_elas)/f_elas ### 计算弹性流动方向张量
beta = 3*mu*dp/sig_eq ### 计算硬化参数
new_sig = sig_elas-beta*s ### 更新应力张量
return as_long_vector(new_sig), as_long_vector(n_elas), beta, dp
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*fe.inner(N_elas, e)*N_elas-2*mu*beta*fe.dev(e)
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = fe.dx(metadata=metadata)
# TODO
a_Newton = fe.inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -fe.inner(as_3D_tensor(sig), eps(u_))*dxm + 1e-10*fe.dot(normal, u_)*ds(2)
def local_project(v, V, u=None):
dv = fe.TrialFunction(V)
v_ = fe.TestFunction(V)
a_proj = fe.inner(dv, v_)*dxm
b_proj = fe.inner(v, v_)*dxm
solver = fe.LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = fe.Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
def global_project(v, V, u=None):
'''
Not working.
AssertionError: Mismatch of quadrature points!
'''
if u is None:
return fe.project(v, V)
else:
u.assign(fe.project(v, V))
P0 = fe.FunctionSpace(mesh, "DG", 0)
p_avg = fe.Function(P0, name="Plastic strain")
V = fe.FunctionSpace(mesh, 'CG', 1)
T_crt = fe.interpolate(fe.Constant(ambient_T), V)
T_pre = fe.interpolate(fe.Constant(ambient_T), V)
T_old = fe.interpolate(fe.Constant(ambient_T), V)
T_old.rename('T', 'T')
v = fe.TestFunction(V)
# If theta = 0., we recover implicit Eulear; if theta = 1., we recover explicit Euler; theta = 0.5 seems to be a good choice.
theta = 0.5
T_rhs = theta*T_pre + (1 - theta)*T_crt
T_bcs = [fe.DirichletBC(V, fe.Constant(ambient_T), bottom)]
class LaserExpression(fe.UserExpression):
def __init__(self, t):
super(LaserExpression, self).__init__()
self.t = t
def eval(self, values, x):
t = self.t
values[0] = 2*P*eta/(np.pi*r**2) * np.exp(-2*((x[0] - x0 - vel*t)**2 + (x[1] - y0)**2) / r**2)
def value_shape(self):
return ()
q_laser = LaserExpression(None)
q_convection = h * (T_rhs - ambient_T)
q_top = q_convection + q_laser
q_surr = q_convection
res_T = rho*Cp/dt*(T_crt - T_pre) * v * fe.dx + k * fe.dot(fe.grad(T_rhs), fe.grad(v)) * fe.dx \
- q_top * v * ds(2) - q_surr * v * ds(3)
files_vtk = glob.glob('data/vtk/pbf/sols' + f"/*")
files_xdmf = glob.glob('data/xdmf/pbf/' + f"/*")
for f in files_vtk + files_xdmf:
os.remove(f)
T_vtk_file = fe.File(f'data/vtk/pbf/sols/T.pvd')
T_vtk_file << T_old
u_vtk_file = fe.File(f'data/vtk/pbf/sols/u.pvd')
u_vtk_file << u
p_vtk_file = fe.File(f'data/vtk/pbf/sols/p.pvd')
p_vtk_file << p_avg
file_results = fe.XDMFFile('data/xdmf/pbf/u.xdmf')
file_results.parameters["functions_share_mesh"] = True
file_results.write(T_old, 0)
file_results.write(u, 0)
file_results.write(p_avg, 0)
plastic_inverval = 10
for i in range(len(ts) - 1):
# for i in range(100):
print(f"step {i + 1}, time = {ts[i + 1]}")
q_laser.t = theta*ts[i] + (1 - theta)*ts[i + 1]
solver_parameters = {'newton_solver': {'maximum_iterations': 20, 'linear_solver': 'mumps'}}
fe.solve(res_T == 0, T_crt, T_bcs, solver_parameters=solver_parameters)
T_pre.assign(T_crt)
print(f"min T = {np.min(np.array(T_pre.vector()))}")
print(f"max T = {np.max(np.array(T_pre.vector()))}")
if (i + 1) % plastic_inverval == 0:
T_crt_array = np.array(T_crt.vector())
T_crt_array = np.where(T_crt_array < args.T_ambient, args.T_ambient, T_crt_array)
T_crt_array = np.where(T_crt_array > args.T_melt, args.T_melt, T_crt_array)
T_old_array = np.array(T_old.vector())
T_old_array = np.where(T_old_array < args.T_ambient, args.T_ambient, T_old_array)
T_old_array = np.where(T_old_array > args.T_melt, args.T_melt, T_old_array)
dT = fe.Function(V)
dT.vector()[:] = T_crt_array - T_old_array
Nitermax, tol = 200, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
A, Res = fe.assemble_system(a_Newton, res, u_bcs)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(fe.Constant((0., 0., 0.)))
niter = 0
nRes = 1.
tol = 1e-8
# while nRes/nRes0 > tol and niter < Nitermax:
# while niter < Nitermax:
while nRes > tol:
print("\n")
fe.solve(A, du.vector(), Res, "mumps")
Du.assign(Du+du)
deps = eps(Du)
print(f"du norm = {np.linalg.norm(np.array(du.vector()))}")
sig_, n_elas_, beta_, dp_ = proj_sig(deps, dT, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
print(f"sig dof = {np.array(sig.vector()).shape}")
print(f"sig norm = {np.linalg.norm(np.array(sig.vector()))}")
print(f"n_elas norm = {np.linalg.norm(np.array(n_elas.vector()))}")
print(f"beta norm = {np.linalg.norm(np.array(beta.vector()))}")
A, Res = fe.assemble_system(a_Newton, res, u_bcs)
nRes = Res.norm("l2")
print(f"Residual: {nRes}")
niter += 1
u.assign(u + Du)
sig_old.assign(sig)
p.assign(p + local_project(dp_, W0))
p_avg.assign(fe.project(p, P0))
T_old.assign(T_crt)
T_vtk_file << T_old
u_vtk_file << u
p_vtk_file << p_avg
file_results.write(T_old, i)
file_results.write(u, i)
file_results.write(p_avg, i)
if __name__ == '__main__':
simulation()
However, the following error occurs in the runtime：
Moving new file over differing existing file:
src: /home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16/error.log.57b0160c06294e748ca68aac62a9c044
dst: /home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16/error.log
backup: /home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16/error.log.old
Backup file exists, overwriting.
------------------- Start compiler output ------------------------
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp: In member function 'virtual void ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const':
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp:124:29: error: 'FE17_C0_Q4' was not declared in this scope
124 | TF0[iq] = fw0 * FE17_C0_Q4[0][iq][iq];
| ^~~~~~~~~~
/tmp/tmpx7xoc3ez/ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16.cpp:125:34: error: 'FE17_C0_Q4' was not declared in this scope
125 | BF0[iq][iq] += TF0[iq] * FE17_C0_Q4[0][iq][iq];
| ^~~~~~~~~~
------------------- End compiler output ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16
---------------------------------------------------------------------------
DijitsoError Traceback (most recent call last)
Cell In[8], line 320
315 file_results.write(p_avg, i)
319 if __name__ == '__main__':
--> 320 simulation()
File ~/plasticity-main/src/utils.py:9, in walltime.<locals>.wrapper(*list_args, **keyword_wargs)
7 def wrapper(*list_args, **keyword_wargs):
8 start_time = time.time()
----> 9 func(*list_args, **keyword_wargs)
10 end_time = time.time()
11 time_elapsed = end_time - start_time
Cell In[8], line 287, in simulation()
283 print(f"du norm = {np.linalg.norm(np.array(du.vector()))}")
285 sig_, n_elas_, beta_, dp_ = proj_sig(deps, dT, sig_old, p)
--> 287 local_project(sig_, W, sig)
288 local_project(n_elas_, W, n_elas)
289 local_project(beta_, W0, beta)
Cell In[8], line 164, in simulation.<locals>.local_project(v, V, u)
162 a_proj = fe.inner(dv, v_)*dxm
163 b_proj = fe.inner(v, v_)*dxm
--> 164 solver = fe.LocalSolver(a_proj, b_proj)
165 solver.factorize()
166 if u is None:
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/solvers.py:39, in LocalSolver.__init__(self, a, L, solver_type)
36 self.L_ufl = L
38 # Wrap as DOLFIN forms
---> 39 a = Form(a)
40 if L is None:
41 # Initialize C++ base class
42 cpp.fem.LocalSolver.__init__(self, a, solver_type)
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/form.py:43, in Form.__init__(self, form, **kwargs)
40 if form_compiler_parameters is None:
41 form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
---> 43 ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
44 mpi_comm=mesh.mpi_comm())
45 ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
47 function_spaces = kwargs.get("function_spaces")
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py:47, in mpi_jit_decorator.<locals>.mpi_jit(*args, **kwargs)
45 # Just call JIT compiler when running in serial
46 if MPI.size(mpi_comm) == 1:
---> 47 return local_jit(*args, **kwargs)
49 # Default status (0 == ok, 1 == fail)
50 status = 0
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/jit/jit.py:97, in ffc_jit(ufl_form, form_compiler_parameters)
95 p.update(dict(parameters["form_compiler"]))
96 p.update(form_compiler_parameters or {})
---> 97 return ffc.jit(ufl_form, parameters=p)
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py:217, in jit(ufl_object, parameters, indirect)
214 kind, module_name = compute_jit_prefix(ufl_object, parameters)
216 # Inspect cache and generate+build if necessary
--> 217 module = jit_build(ufl_object, module_name, parameters)
219 # Raise exception on failure to build or import module
220 if module is None:
221 # TODO: To communicate directory name here, need dijitso params to call
222 #fail_dir = dijitso.cache.create_fail_dir_path(signature, dijitso_cache_params)
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/ffc/jitcompiler.py:130, in jit_build(ufl_object, module_name, parameters)
123 params = dijitso.validate_params({
124 "cache": cache_params,
125 "build": build_params,
126 "generator": parameters, # ffc parameters, just passed on to jit_generate
127 })
129 # Carry out jit compilation, calling jit_generate only if needed
--> 130 module, signature = dijitso.jit(jitable=ufl_object,
131 name=module_name,
132 params=params,
133 generate=jit_generate)
134 return module
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dijitso/jit.py:216, in jit(jitable, name, params, generate, send, receive, wait)
212 lib = load_library(signature, cache_params)
214 if err_info:
215 # TODO: Parse output to find error(s) for better error messages
--> 216 raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
217 % err_info['fail_dir'], err_info)
219 # Return built library and its signature
220 return lib, signature
DijitsoError: Dijitso JIT compilation failed, see '/home/ljr/plasticity-main/src/jitfailure-ffc_form_a4835aa8bc5992adce9bec698822651d96d78c16' for details