Hello everyone,
I am trying to modify yield condition by tweaking @bleyerj 's tutorial Elasto-plastic analysis of a 2D von Mises material.
Given von-mises yield condition :
I want to change it to :
where :-
I have tried to do so in the following way (MWE) :-
from dolfin import *
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
import mshr
import numpy as np
domain = mshr.Rectangle(Point(0,0), Point(1.,1.))
mesh = mshr.generate_mesh(domain, 30)
plot(mesh);
E = 10e3
nu = 0.33
sig0C = 60.
sig0T = 45.
m = sig0C/sig0T
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
Et = E/100.
H = E*Et/(E-Et)
deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
def eps(v):
e = sym(grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sk = sig_elas - s
i111 = tr(sk)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = (1./2.*m)*((m - 1)*i111+ (m+1)*sig_eq) - sig0T - H*old_p
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_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1., 1e-8)
boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
bottom().mark(boundaries, 1)
top().mark(boundaries, 2)
loading = Expression("t", t=0, degree=1)
bc_bottom = DirichletBC(V.sub(1), Constant((0.)), boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
bc_u = [bc_bottom]
def F_ext(v):
return loading*dot(n, v)*ds(2)
a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
Nitermax, tol = 200, 1e-8
Nincr = 1000
load_steps = np.linspace(0, 3000., Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
loading.t = t
A, Res = assemble_system(a_Newton, res, bc_u)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "umfpack")
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, Res = assemble_system(a_Newton, res, bc_u)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
k = p_avg.vector().max()
if k >= 0.02:
break
Initially, it converges nicely but as soon as it reaches the yield condition (i.e. plastic strain will become non-zero in next iteration), it diverges.
I read @nate’s post on
and solution to my problem according to above might be to try easing relaxation parameter.
Therefore, I also tried custom newton solver for the same but ran into Arity mismatch :-
What I tried :-
# Define a Nonlinear problem to assemble the residual and Jacobian
class Problem(NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
# Use a Newton solver with custom damping parameter
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")
self.linear_solver().set_from_options()
def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem, iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
info("Newton damping parameter: %.3e" % theta)
x.axpy(-theta, dx)
F = inner(eps(v), sigma_tang(eps(u_)))*dxm + inner(eps(u_), as_3D_tensor(sig))*dxm - F_ext(u_)
J = derivative(F, u)
problem = Problem(J, F)
solver = CustomSolver()
solver.solve(problem, u.vector())
Error :-
---------------------------------------------------------------------------
ArityMismatch Traceback (most recent call last)
/tmp/ipykernel_720/716830184.py in <module>
26 problem = Problem(J, F)
27 solver = CustomSolver()
---> 28 solver.solve(problem, u.vector())
/tmp/ipykernel_720/2688388353.py in F(self, b, x)
7
8 def F(self, b, x):
----> 9 assemble(self.L, tensor=b)
10
11 def J(self, A, x):
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
200 dolfin_form = form
201 else:
--> 202 dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
203
204 # Create tensor
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
58 return form
59 elif isinstance(form, ufl.Form):
---> 60 return Form(form,
61 form_compiler_parameters=form_compiler_parameters,
62 function_spaces=function_spaces)
/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py in __init__(self, form, **kwargs)
41 form_compiler_parameters = {"external_include_dirs": dolfin_pc["include_dirs"]}
42
---> 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])
/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
48 # Just call JIT compiler when running in serial
49 if MPI.size(mpi_comm) == 1:
---> 50 return local_jit(*args, **kwargs)
51
52 # Default status (0 == ok, 1 == fail)
/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
98 p.update(dict(parameters["form_compiler"]))
99 p.update(form_compiler_parameters or {})
--> 100 return ffc.jit(ufl_form, parameters=p)
101
102
/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
215
216 # Inspect cache and generate+build if necessary
--> 217 module = jit_build(ufl_object, module_name, parameters)
218
219 # Raise exception on failure to build or import module
/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
128
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,
/usr/lib/python3/dist-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
163 elif generate:
164 # 1) Generate source code
--> 165 header, source, dependencies = generate(jitable, name, signature, params["generator"])
166 # Ensure we got unicode from generate
167 header = as_unicode(header)
/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
63 compile_object = compile_coordinate_mapping
64
---> 65 code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
66 prefix=module_name, parameters=parameters, jit=True)
67
/usr/lib/python3/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
140 prefix="Form", parameters=None, jit=False):
141 """This function generates UFC code for a given UFL form or list of UFL forms."""
--> 142 return compile_ufl_objects(forms, "form", object_names,
143 prefix, parameters, jit)
144
/usr/lib/python3/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
183 # Stage 1: analysis
184 cpu_time = time()
--> 185 analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
186 _print_timing(1, time() - cpu_time)
187
/usr/lib/python3/dist-packages/ffc/analysis.py in analyze_ufl_objects(ufl_objects, kind, parameters)
87
88 # Analyze forms
---> 89 form_datas = tuple(_analyze_form(form, parameters)
90 for form in forms)
91
/usr/lib/python3/dist-packages/ffc/analysis.py in <genexpr>(.0)
87
88 # Analyze forms
---> 89 form_datas = tuple(_analyze_form(form, parameters)
90 for form in forms)
91
/usr/lib/python3/dist-packages/ffc/analysis.py in _analyze_form(form, parameters)
181 error("Failed to import tsfc.ufl_utils.compute_form_data when asked "
182 "for tsfc representation.")
--> 183 form_data = tsfc_compute_form_data(form)
184 elif r == "quadrature":
185 # quadrature representation
~/.local/lib/python3.10/site-packages/tsfc/ufl_utils.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_restrictions, do_estimate_degrees)
47 other form compilers based on TSFC.
48 """
---> 49 fd = ufl_compute_form_data(
50 form,
51 do_apply_function_pullbacks=do_apply_function_pullbacks,
/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
405 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
406
--> 407 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
408
409 # TODO: This member is used by unit tests, change the tests to
/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176 for itg in form.integrals():
--> 177 check_integrand_arity(itg.integrand(), arguments, complex_mode)
/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in check_integrand_arity(expr, arguments, complex_mode)
157 key=lambda x: (x.number(), x.part())))
158 rules = ArityChecker(arguments)
--> 159 arg_tuples = map_expr_dag(rules, expr, compress=False)
160 args = tuple(a[0] for a in arg_tuples)
161 if args != arguments:
/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress, vcache, rcache)
34 Return the result of the final function call.
35 """
---> 36 result, = map_expr_dags(function, [expression], compress=compress,
37 vcache=vcache,
38 rcache=rcache)
/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress, vcache, rcache)
97 r = handlers[v._ufl_typecode_](v)
98 else:
---> 99 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
100
101 # Optionally check if r is in rcache, a memory optimization
/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in sum(self, o, a, b)
46 def sum(self, o, a, b):
47 if a != b:
---> 48 raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
49 return a
50
ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('v_0', 'v_1').
Kindly suggest some pointers on how to go about the same. Any thoughts are most welcome.
Thanks.