Dear Fenics team
Thank you for your effort to make Fenics more useful everyday.
I have an electroelasto (axisymmetric) problem, which seems quite simple.
Nan error keeps on. I tried to figure out what went wrong but found no
clue.
Please take a look at my script.
Thank you.
Deok-Kee Choi, Korea.
=================== list ====================================
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 = VectorFunctionSpace(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_vector,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"]="newton" # "snes"
#prm["newton_solver"]["absolute_tolerance"]=1e-8
#prm["newton_solver"]["relative_tolerance"]=1e-7
prm["newton_solver"]["linear_solver"]="cg" # "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()
********************** error message ***************************************
Solving linear system of size 123 x 123 (PETSc Krylov solver).
Solving nonlinear variational problem.
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Newton iteration 0: r (abs) = 1.414e+01 (tol = 1.000e-10) r (rel) =
1.000e+00 (tol = 1.000e-09)
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Applying boundary conditions to linear system.
Solving linear system of size 961 x 961 (PETSc Krylov solver).
PETSc Krylov solver starting to solve 961 x 961 system.
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 linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 1000 iterations (PETSc reason
DIVERGED_ITS, residual norm ||r|| = 1.772912e+01).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: unknown
***
*** DOLFIN version: 1.4.0
*** Git changeset: unknown
***
-------------------------------------------------------------------------
fenics@fenicsproject.org
_______________________________________________
fenics mailing list
fenics@fenicsproject.org
http://fenicsproject.org/mailman/listinfo/fenics
