How to use the solution from the last step as initial guess of a nonlinear variation problem?

Question

I try to solve a examplary nonlinear solid mechanics problem with nonuniform body force. The solution does not converge if a zero initial guess is provided after running some steps. I think using the solution from the last step as an initial guess would be helpful for this problem. Basically, I want to provide a non-zero initial guess at the beginning of each Newton-Raphson step, not only just the first step (which can be found here). Can anyone suggest how to do implement this?

My python code is attached as follows:

from dolfin import *
import os, shutil, math, sys, sympy
import numpy as np

# Create mesh and define function space
mesh = Mesh('Meshdata.xml')
#plot(mesh,interactive=True,axes=True)

savedir = "result/"
if os.path.isdir(savedir):
    shutil.rmtree(savedir)  

d = mesh.topology().dim()
e1 = [Constant([1.,0.]),Constant((1.,0.,0.))][d-2]
e2 = [Constant([0.,1.]),Constant((0.,1.,0.))][d-2]

# Mark boundary subdomians
class top(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[1], 0.256)

top = top()
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
top.mark(boundaries, 1)
ds = Measure("ds")[boundaries]

# Define functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

# Define Dirichlet boundary 
load_min = 0. # load multiplier min value
load_max = 1. # load multiplier max value
load_steps = 101 # number of time steps

utop = Expression(("0", "-0.01*t"), t = 0.0)
bcs = DirichletBC(V, utop, top)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = E/(2*(1 + nu)), E*nu/((1 + nu)*(1 - 2*nu))

# Strain and stress
def eps(v):
    return sym(grad(v))

def sigma(v):
    return 2.0*mu*(eps(v)) + lmbda*tr(eps(v))*Identity(d)

x = SpatialCoordinate(mesh)

def B(u,x):
    return Expression(("0", "pow(1/(x[1]+u1),4)), u1=u.sub(1)")

T = Constant((0.0, 0.0))  # Traction force on the boundary1

# Stored strain energy density
psi = 0.5*inner(sigma(u),eps(u))

# Total potential energy
Pi = psi*dx - dot(B(u,x), u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
Pi_u = derivative(Pi, u, v)

# Compute Jacobian of F
Pi_u_u = derivative(Pi_u, u, du)

# Variational problem and solver
problem_u = NonlinearVariationalProblem(Pi_u, u, bcs, Pi_u_u)
solver_u = NonlinearVariationalSolver(problem_u)

# Solver parameters
prm = solver_u.parameters
info(prm, True)
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-6
prm['newton_solver']['maximum_iterations'] = 2000
prm['newton_solver']['relaxation_parameter'] = 0.5
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'hypre_euclid' # hypre_euclid is the closest alternative of ilu
#prm['newton_solver']['preconditioner'] = 'ilu' # ilu does not work in parallel

load_multipliers = np.linspace(load_min,load_max,load_steps)
forces = np.zeros((len(load_multipliers),2))
file_u = File(savedir+"/u.pvd");
file_s = File(savedir+"/stress.pvd");

for (i_t,t) in enumerate(load_multipliers):
    utop.t = t
    solver_u.solve()

    W = TensorFunctionSpace(mesh, "Lagrange", 2)
    sigma_w = project(2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(d), W)

    file_u << u;
    file_s << sigma_w

    forces[i_t] = np.array([t*0.01,assemble(inner(sigma(u)*e2,e2)*ds(1))])
    np.savetxt(savedir+"/force.txt",forces)

The setting of non-zero guess is

parameters['krylov_solver']['nonzero_initial_guess'] = True

But I have no clear idea of how to use the converged solution from the last step as the initial guess of the next Newton-Raphson step. How to provide a nonzero initial guess for the solver?
I would appreciate any suggestions.

Thanks.

Answer 1

I had a similar question.

http://fenicsproject.org/qa/9536/how-to-set-initial-guess-for-nonlinear-variational-problem

How do you save your previous solution?

Comments

Thank you, Aldenpack. I noticed your question but mine is a little bit different. In my case, I need to provide the numerical solution from last step as an initial guess instead of an analytical function. I'm not quite familiar with how to implement this using python. Do you know how to do that?

---

My idea is to save the solution as a temporary array and feed it to the initial guess in the next step in my current code.

---

u is initialized as a function with zero as its coefficients. Perhaps you could call the following in your loop? (not only before the loop but in it as well)

problem_u = NonlinearVariationalProblem(Pi_u, u, bcs, Pi_u_u)
solver_u = NonlinearVariationalSolver(problem_u)

When you run the solver it returns a modified u with new coefficients. You should be able to just restate the problem with the modified u. Perhaps you just have to restate the problem so the solver knows that u has changed?

Answer 2

Here is some small example code. See also this question.

from dolfin import *
parameters["krylov_solver"]["nonzero_initial_guess"] = True
# parameters["krylov_solver"]["monitor_convergence"] = True

mesh = UnitSquareMesh(256, 256)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("sin(x[0]*12) - x[1]")
a = u*v*dx  + inner(grad(u), grad(v))*dx
L = f*v*dx
A = assemble(a)
b = assemble(L)
solver = KrylovSolver("cg", "ilu")
solver.set_operator(A)
sol = Function(V)

# first run: `sol` (used as initial guess) equals zero everywhere
t0 = time()
solver.solve(sol.vector(), b)
print "Time with zero as initial guess: ", time() - t0

# second run: `sol` now contains the solution that we use as initial guess
t1 = time()
solver.solve(sol.vector(), b)
print "Time with answer as initial guess: ", time() - t1

This will give something like:

Solving linear system of size 66049 x 66049 (PETSc Krylov solver).
Time with zero as initial guess: 2.52104496956
Solving linear system of size 66049 x 66049 (PETSc Krylov solver).
Time with answer as initial guess: 0.0516879558563

I noticed that

parameters["krylov_solver"]["nonzero_initial_guess"] = True

must be given before the solver is initialized

solver = KrylovSolver("cg", "ilu")

or it won't work

Answer 3

This functionality is not supported by NonlinearVariationalSolver.

The linear solve within each Newton iteration is for the direction to the next iterate. There's no reason to expect that the new direction is close to the previous one - unless you made a much too small step in the last iteration.

If you want to try it though, you can implement a Newton solver in a few lines of Python code. You can also try switching to linesearch-Newton (use snes as the nonlinear solver) or use Picard iterations to expand the basin of convergence.