solution continuation not working

Issue #57 resolved
Li Dong created an issue

Hi dolfin-adjoint developers,

Thanks for your efforts to develop this awesome tool. Recently I met a problem when I was trying to implement a continuation strategy, i.e. feed the solution from last iteration to the next iteration. I am aware that this has been demonstrated in one of your examples (stokes_topology). But it somehow failed on me when I was testing a dumb case.

For my dumb test case, I took the poisson-mother.py example and did a slight modification. Now it looks like below. Basically, what I did is that, instead of letting it run straight for 20 iterations and end up with the convergence f_val=1.99e-7, I set a 20-iteration for loop and feed each solution as initial guess to the control variable in each for loop. But this test case got stuck at f_val=3.64e-7 from only the second for loop and l-BFGS kept throwing "REL_REDUCTION_OF_F_<=_FACTR*EPSMCH " in all later for loops.

By the way, I'm using 2016.1.0 version for both fenics and dolfin-adjoint.

from dolfin import *
from dolfin_adjoint import *


set_log_level(ERROR)

# Create mesh, refined in the center
n = 64
mesh = UnitSquareMesh(n, n)

# Define discrete function spaces and funcions
V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)

f = interpolate(Expression("0.0"), W, name='Control')

def forward(f):
  u = Function(V, name='State')
  v = TestFunction(V)

  # Define and solve the Poisson equation to generate the dolfin-adjoint annotation
  F = (inner(grad(u), grad(v)) - f*v)*dx
  bc = DirichletBC(V, 0, "on_boundary")
  solve(F == 0, u, bc)
  return u

# Define functional of interest and the reduced functional
x = SpatialCoordinate(mesh)
d = 1/(2*pi**2)*sin(pi*x[0])*sin(pi*x[1])   # the desired temperature profile

alpha = Constant(1e-6)

for it in range(20):
  adj_reset()

  u = forward(f)

  J = Functional((0.5*inner(u-d, u-d))*dx + alpha/2*f**2*dx)
  control = Control(f)
  rf = ReducedFunctional(J, control)
  f_opt = minimize(rf, bounds=(0.0, 0.8), tol=1.0e-10, options={"gtol": 1.0e-10, "factr": 0.0, "maxiter":1})
  f.assign(f_opt)

Comments (4)

  1. Li Dong reporter

    Somehow when I changed the convergence criteria of the solver like below, everything works...

    f_opt = minimize(rf, bounds=(0.0, 0.8), tol=1.0e-10, options={"ftol": 1e-20,"gtol": 1.0e-15, "factr": 0.0, "maxiter":2})
    
  2. Log in to comment