How to set up initial guess for nonlinear variational problem.

Question

I have a complicated nonlinear system that may have different solutions depending on my initial guess. I want to know how I can use the NonlinearVariationalProblem with an initial guess.

The FEniCS tutorial on page 73 has an example of how to implement a nonlinear variational problem. I have copied it below. The tutorial says that u_ is the most recently computed solution. I'm confused, is this a solution that the solver will find and then reuse or do I need to give it some initial guess? If I don't give it an initial guess what guess will it use?

du = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V) # the most recently computed solution
F = inner(q(u_)nabla_grad(u_), nabla_grad(v))dx
J = derivative(F, u_, du) # Gateaux derivative in dir. of du

problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver = NonlinearVariationalSolver(problem)
...
solver.solve()

How could I make an initial guess so that the system converges easier? (I have an idea more or less of what the solution should look like.) Would I insert some initial solution for u_?

u_ = Function(V)
u_.interpolate(Expression*)

Any input would be very helpful.

Answer 1

Hi, with the above u_ = Function(V) gives you a zero initial guess (Function(V) has a vector of zeros as expansion coefficients). You are right about getting a non-zero initial guess. For completeness here is a modified nonlinear Poisson demo that shows this

from dolfin import * 

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
# k = 0 - zero init guess. See what happens when k is different. k = -1 yields divergence 
u = interpolate(Expression('std::pow(x[0], k)*sin(k*x[1])', k=3), V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})