Hi
I have to solve a variational 2D problem with a mixed function space as:
U=VectorFunctionSpace(mesh,"P",1,2)
Unidfs=FunctionSpace(mesh,"P",1)
Matr=TensorFunctionSpace(mesh, "P", 1)
d=Function(Unidfs)
Trial_d=TrialFunction(Unidfs)
Test_d=TestFunction(Unidfs)
u_buid= VectorElement("P", mesh.ufl_cell(), 1,2)
Fv_build=TensorElement('P',mesh.ufl_cell(),1,(2,2))
el_build=[u_buid,Fv_build]
elment_build=MixedElement(el_build)
Space_vec=FunctionSpace(mesh, elment_build)
state_var=Function(Space_vec)
state_var_Trial=TrialFunction(Space_vec)
state_var_Test=TestFunction(Space_vec)
#############Material's parameters############################
# ######elastic
mu=2
#kelvin springs
beta_vec=1
#viscous
tau_vec=.8
###############damage
p_inf=2
p_vec=2
G_frac=5
zita=1
Du_damage=2
b_int=1#np.real((Du_damage/2)/((2/sqrt(1-zita))*np.arcsinh(sqrt(1-zita)/sqrt(zita))))
c_0_dam=c_0formula(zita)
kappa1=.1
print('mu_inf=',str(mu),'alfa=','beta=',str(beta_vec),'tau=',tau_vec)
print('p_inf=',p_inf,'p=',p_vec,'G_frac=',G_frac,'zita=',zita,'D_u=',Du_damage,'b=',b_int,'c0=',c_0_dam)
#####################evolution parameter################################
tin=0
speed_test=1
massima_deformazione=6
umax=massima_deformazione*lpx
tfin=70#(lambda_lim+5)*lpx
dt=1
t=tin
uL=0
def displacement (dt,predispl,t):
dimp=speed_test*dt+predispl
dpre=interpolate(Constant(0),Unidfs)
####################definitions########################
def F(u): return grad(u) + Identity(2)
def C(u):return (F(u).T)*F(u)
def degrad_func_inf(d) : return (1-d)**p_inf
def degrad_func_i(d,p_i):return (1-d)**p_i
def alpha_dam(d): return zita*d+((1-zita)*(d**2))
def gamma_dam(d): return (1/c_0_dam)*((alpha_dam(d)/b_int)+(b_int*dot(grad(d),grad(d))))
def damage_rate_component(d):
if evol_type=='tensile test':
time_div=dt
elif evol_type=='tensile test with focus':
if t<elastic_path:
time_div=mult_t*dt
else:
time_div=dt
return kappa1*(d-dpre)/time_div
def Energy_inf(T): return (mu/2)*(tr(T)+(1/det(T))-3)
def Energy_i(T,Visc): return (mu/2)*(inner(T,inv(Visc.T*Visc))+(det(inv(Visc.T*Visc))/det(T))-3)
def Total_energy(State_variable,d):
u,Fv=split(State_variable)
psi=(Energy_inf(C(u))*degrad_func_inf(d))+(beta_vec*degrad_func_i(d,p_vec)*Energy_i(C(u),Fv))
return psi
def P_inf_el(State_variable):
u,Fv=split(State_variable)
return mu*(F(u)-(4*det(inv(C(u)))*inv(F(u).T)))
def Q_alpha(State_variable,Q_alpha_pre_vec,P_inf_pre,d):
time_exp=-dt*degrad_func_i(d,p_vec)/(2*tau_vec)
Q=((np.e**(2*time_exp))*Q_alpha_pre_vec)+((np.e**time_exp)*beta_vec*(P_inf_el(State_variable)-P_inf_pre))
return Q
def energy_desity(State_variable,d): return Total_energy(State_variable,d)+(G_frac*gamma_dam(d))+damage_rate_component(d)
##################### setup damage problem#########################ààà
d_min=interpolate(Constant(0),Unidfs)
d_max=interpolate(Constant(1-DOLFIN_EPS),Unidfs)
#boundary conditions Dirichlet (5 lines), if you want use the Newmann b.c., comment the next 5 lines
bcd_dleft=DirichletBC(Unidfs,0,boundary_left)
bc_dright=DirichletBC(Unidfs,0,boundary_right)
bcd=[bcd_dleft,bc_dright]
for bc in bcd:
bc.apply(d_max.vector())
class DamageProblem(OptimisationProblem):
def f(self,x):
d.vector()[:]=x
return assemble(energy_desity(state_var,d)*dx)
def F(self,b,x):
d.vector()[:]=x
assemble(derivative(energy_desity(state_var,d)*dx,d,Test_d),b)
def J(self,A,x):
d.vector()[:]=x
assemble(derivative(derivative(energy_desity(state_var,d)*dx,d,Test_d),d,Trial_d),A)
solver_d=PETScTAOSolver()
solver_d.parameters.update({'maximum_iterations':400,"gradient_relative_tol" : 1.0e-2,'report':True})
#solver_d.parameters.update({'method':'tron','linear_solver':'umfpack','line_search':'gpcg','maximum_iterations':400})
##################### before start evolution and iterations##########
tol=1e-7
maxite=100
Q_alpha_pre_vec=project(Identity(2),Matr)
P_inf_pre=project(Identity(2)-Identity(2),Matr)
############ start evolutions
while t<tfin:
print(name_dir)
##############################update b.c. and evolution test's parameters####################
err=1
ite=0
bc_left=DirichletBC(Space_vec.sub(0), Constant((0,0)), boundary_left)
bc_right=DirichletBC(Space_vec.sub(0),Constant((uL,0)),boundary_right)
bc=[bc_left,bc_right]
while err>tol and ite<maxite:
print('iter=',ite,'t=',t)
############ define and solve the newton-raphson solver (tangent method) for u at d fixed###############
Energy=energy_desity(state_var,d)*dx
cons_force=derivative(Energy,state_var,state_var_Test)
u_test,Fv_test=split(state_var_Test)
non_cons=inner(Q_alpha(state_var,Q_alpha_pre_vec,P_inf_pre,d),Fv_test)*dx
var1=cons_force-non_cons
var2=derivative(var1,state_var,state_var_Trial)
problem = NonlinearVariationalProblem(var1,state_var,bc,var2)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['maximum_iterations'] = 20
solver.parameters['newton_solver']['linear_solver'] ='umfpack'# 'default'#'cg'#'gmres'#'mumps'#
#solver.parameters.update(snes_solver_parameters)
solver.solve()
################# solve the bounded minimization problem for d at u fixed###############
solver_d.solve(DamageProblem(),d.vector(),d_min.vector(),d_max.vector())
################update error and iter and upre and dpre##################
deltad=project(d-dpre,Unidfs)
err=norm(deltad.vector(),'linf',mesh)#+(norm(deltau.vector(),'linf',mesh)/100)
print(err)
ite=ite+1
dpre.assign(d)
d_min.assign(project(d-deltad,Unidfs))
P_inf_pre.assign(project(P_inf_el(state_var),Matr))
Q_alpha_pre_vec.assign(project(Q_alpha(state_var,Q_alpha_pre_vec,P_inf_pre,d),Matr))
excuse me for the long code, but it gives the following error on solver.solve()
UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory
Blockquote
*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /build/dolfin-DneRIZ/dolfin-2019.2.0~git20201207.b495043/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------