How to avoid repeated assembly correctly in dolfin-adjoint?

Issue #63 resolved
Charles Zhang created an issue

Hi there,

I was following the FEniCS tutorial on avoiding-assembly to assemble and keep a matrix, beforing solving it with a rhs in each time iteration, but dolfin-adjoint threw an error:

Exception TypeError: "'NoneType' object is not callable" in <bound method KeyedDict.__del__ of {(Form([Integral(Sum(Sum(Product(IntValue(-1), Product(Indexed
...

A minimal code is attached, but it gives this error if method 1 is used. Could you advice the correct way to assemble a matrix in dolfin-adjoint environment? I didn't rememble this error with dolfin-adjoint 1.5.

from dolfin import *
from dolfin_adjoint import *

nx = 30
mesh = UnitSquareMesh(nx, nx)

BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG  = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)

def main(f, method=1, annotate=False):
    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)

    w0 = Function(W)
    (u0, p0) = split(w0)

    timestep = Constant(1.0e-1)
    t = 0.0
    end = 1.0e0

    a = (dot(u, v) - div(v)*p + timestep*div(u)*q + p*q)*dx
    L = timestep*f*q*dx + p0*q*dx

    if method==1:
        A = assemble(a) # assemble matrix once

    w = Function(W)

    if annotate: adj_start_timestep()

    while t < end - 0.5*float(timestep):

        print 't = ', t

        if method==1:
            b = assemble(L) # assemble rhs
            solve(A, w.vector(), b, annotate=annotate)
        else:
            solve(a==L, w, annotate=annotate)
        w0.assign(w)
        t += float(timestep)
        if annotate: adj_inc_timestep(t, t > end - .5*float(timestep))
    return w

if __name__ == '__main__':

    f = Constant(-1.0)
    w = main(f, method=1, annotate=True)
    # method 2 works fine,
    # while method 1 throws an error:
    # Exception TypeError: "'NoneType' object is not callable"
    u, p = split(w)

    J = Functional(dot(p, p)*dx*dt[FINISH_TIME])
    dJdf = compute_gradient(J, Control(f))
    print 'dJdf = ', float(dJdf)

Comments (2)

  1. Log in to comment