Switching from Constant to TimeConstantParameter causes a filed taylor test

[Gabriel Balaban]

I have slightly modified the checkpointing tutorial so that the constant nu is now a dolfin scalar Function. The values that the Taylor test returned were

Convergence orders for Taylor remainder with gradient information (should all be 2): [1.062702948792387, 1.1094408461147007, 1.2031069364153415, 1.3571291965740102]

Is there a way to fix this?

Cheers,
Gabriel

from dolfin import *
from dolfin_adjoint import *

adj_checkpointing(strategy='multistage', steps=11,
                  snaps_on_disk=2, snaps_in_ram=2, verbose=True)

n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)

ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])")), V)

def main(nu):
  u = Function(ic)
  u_next = Function(V)
  v = TestFunction(V)

  timestep = Constant(0.01)

  F = (inner((u_next - u)/timestep, v)
     + inner(grad(u_next)*u_next, v)
     + nu*inner(grad(u_next), grad(v)))*dx

  bc = DirichletBC(V, (0.0, 0.0), "on_boundary")

  t = 0.0
  end = 0.1
  while (t <= end):
    solve(F == 0, u_next, bc)
    u.assign(u_next)
    t += float(timestep)
    adj_inc_timestep()

  return u

if __name__ == "__main__":
  VV = FunctionSpace(mesh, "CG", 2) #Modified
  nu = interpolate(Constant(0.0001), VV, name = "nu") #Modified
  u = main(nu)

  J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
  dJdnu = compute_gradient(J, TimeConstantParameter("nu"), forget = False) #Modified

  Jnu = assemble(inner(u, u)*dx)

  parameters["adjoint"]["stop_annotating"] = True
  def Jhat(nu):
    u = main(nu)
    return assemble(inner(u, u)*dx)

  conv_rate = taylor_test(Jhat, TimeConstantParameter("nu"), Jnu, dJdnu) #Modified

[Patrick Farrell (pefarrell)]

Hi Gabriel,

Yep. The fix is simple.

The clue is the Taylor remainders. Even the first-order ones aren't first order:
Convergence orders for Taylor remainder without gradient information (should all be 1): [0.5612058296487964, 0.4863431171924579, 0.4194999069935161, 0.4390197925250572]

This always means that your seed value for the Taylor test is too large. Those numbers have nothing to do with dolfin-adjoint, and should always be 1 (even if your model isn't differentiable, say). As discussed on http://dolfin-adjoint.org/documentation/verification.html, the fix is to pass a smaller seed value. Passing seed=1.0e-5 to taylor_test yields:
Convergence orders for Taylor remainder with gradient information (should all be 2): [1.9882784834167722, 1.9941105578851968, 1.9970480600786011, 1.9985222301337167]

This makes sense because the seed value depends on the scaling of the parameter. If the parameter is a ScalarParameter, taylor_test scales the seed from that one number. But when the parameter is a field, it can't (it doesn't know your field happens to be spatially constant). That's why changing ScalarParameter -> TimeConstantParameter made the difference.

Hope this helps!

Patrick

[Gabriel Balaban (gabrib)]

Thanks a lot Patrick! That was very well explained. Changing the seed value does the trick.
Cheers,
  Gabriel