Initial guess in NonLinearVariationalProblem

Hi,

Does anyone know how NonlinearVariationalProblem sets the initial guess for the newton iteration at each load step when the the Dirichlet B.C. is changed?

I can reproduce the first load step in my own implementation of the newton method and the initial guess appears to be from solving the linearized problem with the first initial guess. However, in subsequent changes of Dirichlet B.C. (u_R and u_up), the initial guess doesn't seems to come from the solution of the previous load step. It also appear that the first iteration (of the 2nd load step) in NonlinearVariationalProblem is twice of that from my own implementation. Any idea?

Thanks!

from dolfin import *
import numpy as np

parameters["form_compiler"]["representation"] ="uflacs"
parameters["form_compiler"]["quadrature_degree"]= 2


def Fmat(u):

    d = u.ufl_domain().geometric_dimension()
    I = Identity(d)
    F = I + grad(u)
    return F

#boundary conditions
class Left(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[0]) < tol

class Right(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[0]-10.0) < tol

class Lower(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[1]) < tol

class Upper(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[1]-10.0) < tol

class Bottom(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[2]) < tol

class Top(SubDomain):
      def inside(self, x, on_boundary):
          tol = 1E-14
          return on_boundary and abs(x[2]-1.0) < tol

nx = 1
ny = 1
nz = 1
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10.0, 10.0, 1.0), nx, ny, nz)
V = VectorFunctionSpace(mesh, 'CG', 2)
TF = TensorFunctionSpace(mesh,'CG',1)
Q = FunctionSpace(mesh,'CG',1)
W = MixedFunctionSpace([V,Q])

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left = Left()
right = Right()
bottom = Bottom()
top = Top()
lower = Lower()
upper = Upper()

left.mark(boundaries, 1)
right.mark(boundaries, 2)
lower.mark(boundaries, 3)
upper.mark(boundaries, 4)
bottom.mark(boundaries, 5)
top.mark(boundaries, 6)

c = Constant(0.0)
u_R = Expression('x[0]*(lmbda - 1.0)', lmbda = 1.0)
u_up = Expression('x[1]*(lmbda - 1.0)', lmbda = 1.0)

bcleft0 = DirichletBC(W.sub(0).sub(0), c, boundaries, 1)
bcright0 = DirichletBC(W.sub(0).sub(0), u_R, boundaries, 2)
dbcright0 = DirichletBC(W.sub(0).sub(0), c, boundaries, 2)
bclower1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 3)
bcupper1 = DirichletBC(W.sub(0).sub(1), u_up, boundaries, 4)
dbcupper1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 4)
bcbottom2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 5)

bcs = [bcleft0, bcright0, bclower1, bcupper1, bcbottom2]
bcs_dw = [bcleft0, dbcright0, bclower1, dbcupper1, bcbottom2]

dw = TrialFunction(W)
du,dq = TrialFunctions(W)
v,q = TestFunctions(W)
wtest = TestFunction(W)
w = Function(W)
u,p = split(w)


J = det(Fmat(u))
I1 = tr(Fmat(u).T*Fmat(u))
Wsimple = ((I1**2 - 9.0) - p*(J - 1.0))*dx

F1 = derivative(Wsimple,w,wtest)

Jac = derivative(F1, w, dw)
############################ Using Own Newton Method ############################################
abs_tol = 1.0E-10
rel_tol = 1.0E-9
maxiter = 25
omega = 1.0
sigma_real = []
lmbda = 1.0
dxx = 1.0

w_k = Function(W)
for lmbda_value in range(1, 3):

    it = 0
        print "lmbda value = ", lmbda_value

        u_R.lmbda = 1.0 + lmbda_value * 0.5
        u_up.lmbda = 1.0 + lmbda_value * 0.5


            A0, b0 = assemble_system(Jac, -F1, bcs)
            solve(A0, w.vector(), b0, "mumps")
    resid0 = b0.norm("l2")
    rel_res = b0.norm("l2")/resid0
    res = resid0
    print "Residual: %.3e, Relative residual: %.3e" %(res, rel_res)


    dww = Function(W)

    while (rel_res > rel_tol and res > abs_tol) and it < maxiter:
                    it += 1

                    A, b = assemble_system(Jac, -F1, bcs_dw)
                    solve(A, dww.vector(), b, "mumps")

            B = assemble(F1)
            for bc in bcs_dw:
                    bc.apply(B)

            rel_res = B.norm("l2")/resid0
            res = B.norm("l2")
            w.vector().axpy(1.0, dww.vector())

            print "Residual: %.3e, Relative residual: %.3e" %(res, rel_res)



############################ Using Nonlinearvariationalproblem ############################################
w.vector()[:] = 0
problem = NonlinearVariationalProblem(F1, w, bcs, Jac)
solver = NonlinearVariationalSolver(problem)

displacementfile = File("displacement.pvd")
stressfile = File("stress.pvd")

for lmbda_value in range(1, 3):

        print "lmbda value = ", lmbda_value

        u_R.lmbda = 1.0 + lmbda_value * 0.5
        u_up.lmbda = 1.0 + lmbda_value * 0.5

        solver.solve()