Hi, everyone.
I need to solve a simple electroelastic problem but have NaN errors. The problem is of finding the displacement $u$ of hyperelastic material subjected to electric potential. The unknowns are the displacement $u$ and the potential $\phi$. I don't know what went wrong.
Thank you.
I've got the error message:
1 SNES Function norm nan *** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_FNORM_NAN. Traceback (most recent call last): File "/home/dave/Documents/Research/ElectroElasto/electroelasto_mine_axi_1002.py", line 152, in <module>
solver.solve() 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@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with PETScSNESSolver.
*** Reason: Solver did not converge.
*** Where: This error was encountered inside PETScSNESSolver.cpp.
*** Process: unknown
***
*** DOLFIN version: 1.4.0
*** Git changeset: unknown
*** ------------------------------
=========================================================Soruce==========
from dolfin import *
import numpy as np
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Geometry
ri=10.0
ro=20.0
zb=0.0
zt=1.0
domain = Rectangle(ri,zb,ro,zt)
# Generate and plot mesh
mesh = Mesh(domain,20)
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],ri) and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],ro) and on_boundary
# Initialize sub-domain instances
left = Left()
right=Right()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries,1)
right.mark(boundaries,2)
# Create mesh and define function space
V1 = FunctionSpace(mesh, "Lagrange",2)
V2 = FunctionSpace(mesh,"Lagrange",1)
V= MixedFunctionSpace([V1,V2])
u=TrialFunction(V1)
phi=TrialFunction(V2)
vq=TestFunction(V)
(v,q)=TestFunctions(V)
dup=TrialFunction(V)
up=Function(V)
(u,phi)=split(up)
#v_u1=V1.sub(0)
#v_u2=V1.sub(1)
plot(mesh)
# Define Dirichlet boundary (x = 0 or x = 1)
zero=Constant(0.0)
zero_vector=Constant((0.0,0.0))
volt=Constant(5.0)
bcul = DirichletBC(V1,zero,left)
#bcur = DirichletBC(V1,zero,right)
bcpl = DirichletBC(V2,volt,left)
bcpr = DirichletBC(V2,-volt,right)
bcs = [bcul,bcpl,bcpr]
# Define functions
B = Constant((0.0,0.0)) # Body force per unit volume
T = Constant((0.0,0.0)) # Traction force on the boundary
Q = Constant(0.0) # Traction force on the boundary
r=Expression("x[0]")
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
#F = I + grad(u) # Deformation gradient
F=as_matrix([[1.0+u.dx(0),0.0],[0.0,1.0]])
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
Jf = det(F)
# Elasticity parameters
mu, k = 5.0, 10.0
lmbda = Constant(k-2.0*mu/3.0)
c1, c2 = 10.0,6.0
#EE=-grad(phi)
EE=as_vector([-phi.dx(0),0.0])
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - d) - mu*ln(Jf) + (lmbda/2)*(ln(Jf))**2\
+c1*inner(I,outer(EE,EE))+c2*inner(C,outer(EE,EE))
# Total potential energy
Pi=2.0*np.pi*r*psi*dx
#Pi=psi*dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
Func = derivative(Pi, up, vq)
# Compute Jacobian of F
J = derivative(Func, up, dup)
# Solve variational problem
#solve(Func == 0, up, bcs, J=J,form_compiler_parameters=#ffc_options)
problem=NonlinearVariationalProblem(Func,up,bcs=bcs,J=J,\
form_compiler_parameters=ffc_options)
solver=NonlinearVariationalSolver(problem)
prm=solver.parameters
prm["nonlinear_solver"]="snes" # "snes"
#prm["newton_solver"]["absolute_tolerance"]=1e-8
#prm["newton_solver"]["relative_tolerance"]=1e-7
prm["newton_solver"]["linear_solver"]="gmres" # "cg" "gmres"
prm["newton_solver"]["krylov_solver"]["relative_tolerance"]=1e-6
prm["newton_solver"]["krylov_solver"]["maximum_iterations"]=1000
#prm["newton_solver"]["krylov_solver"]["monitor_convergence"]=True
prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"]=False
#prm["newton_solver"]["krylov_solver"]["gmres"]["restart"]=40
prm["newton_solver"]["preconditioner"]="ilu" # "ilu" "jacobi"
prm["newton_solver"]["krylov_solver"]["preconditioner"]["structure"]\
="same_nonzero_pattern"
#prm["newton_solver"]["krylov_solver"]["preconditioner"]["ilu"]["fill_level"]=0
#solver.parameters["linear_solver"]="cg"
#solver.parameters["preconditioner"]="ilu"
#solver.parameters["krylov_solver"]["absolute_tolerance"]
#solver.parameters["krylov_solver"]["relative_tolerance"]
set_log_level(PROGRESS)
solver.solve()
# Save solution in VTK format
(u,phi)=up.split()
file = File("displacement.pvd")
file << u
file = File("potential.pvd")
file << phi
# Plot and hold solution
plot(u,title="displacement",mode="displacement")
#plot(u.sub(1),title="y-displacement")
plot(phi,title="potential")
interactive()
