I made the changes suggested by you in the pdf but now the equation is solved in 0 iterations as initial guess is (0, 0, 0). When I am trying non-zero guess like (0.01, 0., 0.), the newton solver converges but I am not able to make any useful interpretation from the results.
Can you please tell me where I am going wrong ?
Here is the modified code -
from dolfin import *
from ufl import indices
xL=20.0
yL=20.0
zL=20.0
mesh = BoxMesh(Point(-xL /2 , -yL /2 , -zL ) , Point( xL/2 , yL/2 , 0), 5, 5, 5)
V = VectorFunctionSpace(mesh , 'P' , 1)
T = TensorFunctionSpace(mesh , 'P' , 1)
left = CompiledSubDomain( 'near( x[0] , l ) && on_boundary' , l = -xL/2)
right = CompiledSubDomain( 'near( x[0] , l ) && on_boundary' , l=xL/2)
back = CompiledSubDomain( 'near( x[1] , l ) && on_boundary' , l=-yL/2)
front= CompiledSubDomain( 'near( x[1] , l ) && on_boundary' , l=yL/2)
bottom = CompiledSubDomain( 'near( x[2] , l ) && on_boundary' , l=-zL)
top = CompiledSubDomain( 'x[2] == 0.0 && (x[0]*x[0]) + (x[1]*x[1]) <= r*r ', r=0)
cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
dA = Measure( 'ds' , domain=mesh , subdomain_data= facets, metadata = {"quadrature_degree": 2})
dV = Measure( 'dx' , domain=mesh , subdomain_data= cells, metadata = {"quadrature_degree": 2})
facets.set_all(0)
nanoindent = Expression((' 0.0' , '0.0' , '-Amp' ) , degree = 1, Amp = 1.)
bc = [ DirichletBC(V, ( 0. , 0. , 0. ) , bottom ) ]
f_gr = Constant((0. , 0. ,0. ))
du = TrialFunction(V)
delu = TestFunction(V)
u00 = Function(V)
u0 = Function(V)
u = Function(V)
S=Function(T)
S0=Function(T)
t = 0.0
tend = 3000.0
dt = 100.
init = Expression(( '0.' , '0.' , '0.' ), degree = 0)
u.interpolate(init )
u0.assign(u)
u00.assign(u0)
rho0 = 9000.0*(10**-15)
lambada = 90.0
E1 , E2 = 200.0 , 200.0
mu = 2.0*(10**5)
i,j,k,r = indices(4)
delta = Identity(3)
F = as_tensor(delta[i, j] + u[i].dx(j), [i, j])
eps= as_tensor( 1.0/2.0*(u[i].dx(k) + u[k].dx(i)), [i, k])
eps0= as_tensor( 1.0/2.0*(u0[i].dx(k) + u0[k].dx(i)), [i,k])
eps_dev= as_tensor(eps[i, j] - 1.0/3.0*eps[k, k]*(delta[i, j]), [i, j])
eps0_dev= as_tensor(eps0[i, j] - eps0[k, k]*1.0/3.0*(delta[i, j]), [i, j])
devS0= as_tensor(S0[i, j] - S0[k, k]*1.0/3.0*delta[i, j], [i, j])
S = as_tensor(lambada*eps[ j, j ]*delta[ i, k ] + mu/(E2*dt+mu)*devS0[ i, k] + E1*E2*dt / (E2*dt+mu)*eps_dev[ i, k ] + mu*( E1+E2 ) /( E2*dt+mu)*(eps_dev[i, k] - eps0_dev[i, k]), [i, k])
N = FacetNormal(mesh)
Form = (rho0/dt/dt*(u - 2.*u0 + u00)[i]*delu[i] + F[i, k]*S[r, k]*delu[i].dx(r) - rho0*f_gr[i]*delu[i])*dV - nanoindent[i]*delu[i]*dA(1)
nz = as_tensor([ 0.0 , 0.0 , 1.0])
forceZ = as_tensor(F[i, k]*S[r, k]*nz[r], [i,])
Gain = derivative(Form , u , du )
file_u = File('disp.pvd')
u_max = 0.1
A = 24.5*u_max**2
radius = sqrt(A/pi)
r = radius
top.mark(facets , 1)
time_values =[0.0]
u_max_values =[0.0]
forceZ_values =[0.0]
solver_parameters = {
"newton_solver": {
"linear_solver": "cg",
"preconditioner" : "hypre_amg",
"relative_tolerance" : 1e-2,
"absolute_tolerance" : 1e-5,
"maximum_iterations" : 30
}
}
form_compiler_parameters = {
"cpp_optimize": True,
"representation": "quadrature",
"quadrature_degree": 2
}
while t<tend:
t += dt
if t < 1000.: nanoindent.Amp = 0.002*t/A
if t >= 1000. and t <=1500.: nanoindent.Amp = 2./A
if t > 1500.: nanoindent.Amp = (2. - 0.002*(t - 1500.))/A
if t >=2500.: nanoindent.Amp = 0.
solve(Form == 0 , u , bc , J=Gain , solver_parameters= solver_parameters,form_compiler_parameters= form_compiler_parameters)
file_u << (u,t)
time_values.append(t)
u_max = abs(u((0 ,0 ,0))[2])
u_max_values.append(u_max)
fZ = project(forceZ,V)
fZvalue = abs(fZ((0 ,0 ,0))[2])*A
forceZ_values.append(fZvalue)
S0.assign(project(S ,T))
u00.assign(u0)
u0.assign(u)
Newton-solver on initial trivial guess -
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-05) r (rel) = -nan (tol = 1.000e-02)
Newton solver finished in 0 iterations and 0 linear solver iterations.
PS - I have not refined the mesh as my memory keeps running out when I do it.
