Taylor test for time dependent problem

Issue #81 new
Jacob Downs created an issue

Hello, I am trying to conduct a Taylor test for a simple time dependent problem. The following code appears to work:

from dolfin import *
from dolfin_adjoint import *

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

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

def main(nu):
    u = ic.copy(deepcopy=True)
    u_next = Function(V)
    v = TestFunction(V)

    timestep = Constant(0.005)

    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
    Jnu = 0.0
    end = 0.1
    times = [t]
    adj_start_timestep(time=t)
    while (t <= end):
        solve(F == 0, u_next, bc)
        u.assign(u_next)
        Jnu += assemble(inner(u,u)*dx)
        t += float(timestep)
        times.append(t)
        adj_inc_timestep(time=t, finished=t>end)

    return u,times, Jnu

if __name__ == "__main__":
    nu = Constant(0.0001)

    u,times, Jnu = main(nu)
    Jfull = sum([inner(u,u)*dx*dt[t] for t in times])
    J = Functional(Jfull)

    #J = Functional(inner(u, u)*dx*dt)

    dJdnu = compute_gradient(J, Control(nu))
    print "Gradient is: ", float(dJdnu)

    parameters["adjoint"]["stop_annotating"] = True # stop registering equations

    def Jhat(nu): # the functional as a pure function of nu
        u,times, Jnu = main(nu)
        return Jnu

    conv_rate = taylor_test(Jhat, Control(nu), Jnu, dJdnu)

However, if I change the functional to the following

J = Functional(inner(u, u)*dx*dt)

and change Jnu to this:

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

I do not get the expected rate of convergence. This is a pretty simple integration scheme, but I thought it would work if the timestep was small. How should I be calculating Jnu?

Thank You!

Jacob

Comments (1)

  1. Simon Funke

    Hi Jacob,

    The mistake is that the dt integral implements a trapezoidle rule, whic you implemented a right Riemann sum. In other words, you need to update your Jnu calculation to get the second order convergence again.

    Best wishes,

    Simon

  2. Log in to comment