diff --git a/README b/README
new file mode 100644
index 0000000..5cd8882
--- /dev/null
+++ b/README
@@ -0,0 +1,16 @@
+# Dependencies:
+# Python3
+# fenicsproject (stable or dev)
+# numpy
+
+# fenicsproject run stable
+# fenicsproject run dev
+fenicsproject run quay.io/dolfinadjoint/pyadjoint:latest
+
+# Module files
+cases.py
+optimalflow.py
+
+# Run test script
+run_case.py
+
diff --git a/cases.py b/cases.py
new file mode 100644
index 0000000..72a6678
--- /dev/null
+++ b/cases.py
@@ -0,0 +1,212 @@
+import os
+import os.path
+import csv
+import numpy
+
+from dolfin import *
+from optimalflow import non_negativize
+
+def hdf5write(name, field):
+    mesh = field.function_space().mesh()
+    file = HDF5File(mesh.mpi_comm(), name, "w")
+    file.write(field, "/function", 0)
+    file.close()
+
+def csvwrite(name, values, header=None, debug=True, mode="w"): 
+    if debug:
+        print(header)
+        print(" & ".join("$%.4f$" % v for v in values))
+
+    with open(name, mode) as f:
+        writer = csv.writer(f)
+        writer.writerow(header)     
+        writer.writerow(values)     
+
+def generate_data_bump(mesh, g, D=0.0, r=0.0, noise=0.0, output=None):
+    """Generate data for a relatively smooth test case with a bump moving
+    in the middle of the domain.
+
+    Domain: [0, 1]**2
+    c_0: Bump function 
+    phi = (g, g)
+
+    Let c_1 be defined by one timestep of the convection-diffusion-reaction equation
+
+      c_t + div (phi c) - D div grad c + r c = 0    
+    
+    where c_0 and phi are as given, and D and r if given. Evolve by a
+    timestep tau = 1.0.
+
+    If outputdir is given as a string, write generated data to given
+    directory.
+    """
+    #tau = Constant(1.0) # Time step
+    tau = 1.0 # Time step
+
+    Q = FunctionSpace(mesh, "CG", 1)
+
+    # Define initial c0
+    bump_x = "beta/(2*alpha*tgamma(1./beta))*std::exp(-pow((std::abs(x[0] - mu)/alpha), beta))"  
+    bump_y = "beta/(2*alpha*tgamma(1./beta))*std::exp(-pow((std::abs(x[1] - mu)/alpha), beta))"  
+    bump = Expression("1 + 1./7*%s*%s" % (bump_x, bump_y), degree=3, beta=8.0, alpha=0.2, mu=0.5)
+    c0 = interpolate(bump, Q)
+    
+    # Define initial phi 
+    V = VectorFunctionSpace(mesh, "CG", 1)
+    phi = Expression((str(g), str(g)), degree=1)
+    phi = interpolate(phi, V)
+
+    # Define variational problem for evolution of c
+    c = TrialFunction(Q)
+    d = TestFunction(Q)
+    a = (inner(c, d) + tau*inner(div(c*phi), d))*dx()  
+    if D is not None and abs(D) > 0.0:
+        a += tau*inner(D*grad(c), grad(d))*dx()
+    if r is not None and abs(r) > 0.0:
+        a += tau*r*c*d*dx()
+
+    L = inner(c0, d)*dx()
+    A = assemble(a)
+    b = assemble(L)
+
+    if D is not None and D > 0:
+        bc = DirichletBC(Q, c0, "on_boundary")
+        bc.apply(A, b)
+    
+    # Solve it
+    c1 = Function(Q)
+    solve(A, c1.vector(), b, "mumps")
+
+    # Add noise to the concentrations:
+    if noise: 
+        N = c0.vector().size()
+        r0 = numpy.random.normal(0, noise, N)
+        r1 = numpy.random.normal(0, noise, N)
+        c0.vector()[:] += r0
+        c1.vector()[:] += r1
+
+    # Zero any negative values 
+    c0 = non_negativize(c0)
+    c1 = non_negativize(c1)
+
+    # Optional: Store solutions to .h5 and .pvd format 
+    if output is not None:
+        filename = lambda f: os.path.join(output, "data", f)
+        hdf5write(filename("c0_e.h5"), c0)
+        hdf5write(filename("c1_e.h5"), c1)
+        hdf5write(filename("phi_e.h5"), phi)
+
+        file = File(filename("c_e.pvd"))
+        file << c0
+        file << c1
+        file = File(filename("phi_e.pvd"))
+        file << phi
+
+        csvwrite(filename("info_e.csv"), (float(tau), D, g), header=("tau", "D", "g"))
+
+    return (c0, c1, phi, tau)
+
+def generate_data_bdry(mesh, g, D=0.0, r=0.0, noise=0.0, tau=1.0, output=None):
+    """Generate data for a test case with a concentration appearing on the boundary.
+
+    Domain: Given
+    c_0: 1.0 on the left and lower boundary if d = 2, on_boundary iff d=3
+    phi = (g, 0) if d == 2, (g, 0, 0) iff d == 3
+
+    Let c_1 be defined by one timestep of the convection-diffusion-reaction equation
+
+      c_t + div (phi c) - D div grad c + r c = 0    
+    
+    where c_0 and phi are as given, and D and r if given. Evolve by a
+    timestep tau = 1.0.
+
+    If outputdir is given as a string, write generated data to given
+    directory.
+    """
+    tau = Constant(tau) # Time step
+
+    Q = FunctionSpace(mesh, "CG", 1)
+    dim = mesh.topology().dim()
+    if dim == 2:
+        bc = DirichletBC(Q, 1.0, "near(x[0], 0.0) || near(x[1], 0.0)", "pointwise")
+    else:
+        bc = DirichletBC(Q, 1.0, "on_boundary")
+    c_ = Function(Q)
+    bc.apply(c_.vector())
+    
+    # Define initial phi 
+    V = VectorFunctionSpace(mesh, "CG", 1)
+    if dim == 2:
+        phi = Expression((str(g), "0.0"), degree=1)
+    else:
+        phi = Expression((str(g), "0.0", "0.0"), degree=1)
+    phi = interpolate(phi, V)
+
+    # Define variational problem for evolution of c
+    c = TrialFunction(Q)
+    d = TestFunction(Q)
+    a = (inner(c, d) + tau*inner(div(c*phi), d))*dx()  
+    if D is not None and abs(D) > 0.0:
+        a += tau*inner(D*grad(c), grad(d))*dx()
+    if r is not None and abs(r) > 0.0:
+        a += tau*r*c*d*dx()
+
+    L = inner(c_, d)*dx()
+    A = assemble(a)
+    b = assemble(L)
+
+    if D is not None and D > 0:
+        bc = DirichletBC(Q, c_, "on_boundary")
+        bc.apply(A, b)
+    
+    # Solve it for c0
+    c0 = Function(Q)
+    solve(A, c0.vector(), b, "mumps")
+
+    # Evolve further to c1
+    c_.assign(c0)
+    b = assemble(L)
+    if D is not None and D > 0:
+        bc = DirichletBC(Q, c0, "on_boundary")
+        bc.apply(A, b)
+
+    c1 = Function(Q)
+    solve(A, c1.vector(), b, "mumps")
+    
+    # Add noise to the concentrations:
+    if noise: 
+        N = c0.vector().size()
+        r0 = numpy.random.normal(0, noise, N)
+        r1 = numpy.random.normal(0, noise, N)
+        c0.vector()[:] += r0
+        c1.vector()[:] += r1
+
+    # Zero any negative values 
+    c0 = non_negativize(c0)
+    c1 = non_negativize(c1)
+
+    # Optional: Store solutions to .h5 and .pvd format 
+    if output is not None:
+        filename = lambda f: os.path.join(output, f)
+        hdf5write(filename("c0_e.h5"), c0)
+        hdf5write(filename("c1_e.h5"), c1)
+        delta_c = project(c1 - c0, c0.function_space())
+        hdf5write(filename("delta_c_e.h5"), delta_c)
+        hdf5write(filename("phi_e.h5"), phi)
+
+        file = File(filename("c_e.pvd"))
+        file << c0
+        file << c1
+        file << delta_c
+        file = File(filename("phi_e.pvd"))
+        file << phi
+
+        csvwrite(filename("info_e.csv"), (float(tau), g, D, r), header=("tau", "g", "D", "r"))
+
+    return (c0, c1, phi, tau)
+
+if __name__ == "__main__":
+
+    mesh = UnitSquareMesh(32, 32)
+    (c0, c1, phi, tau) = generate_data_bump(mesh, 0.02, output="tmp-bump")
+    (c0, c1, phi, tau) = generate_data_bdry(mesh, 0.02, output="tmp-bdry")
diff --git a/optimalflow.py b/optimalflow.py
new file mode 100644
index 0000000..556999f
--- /dev/null
+++ b/optimalflow.py
@@ -0,0 +1,160 @@
+import numpy
+import pylab
+from dolfin import *
+
+def non_negativize(I):
+    """Given a scalar Function c, return the field with negative values
+    replaced by zero.
+
+    """
+    c = I.copy(deepcopy=True) 
+    vec = c.vector().get_local()
+    # Extract dof indices where the vector values are negative
+    negs = numpy.where(vec < 0)[0]
+    # Set these values to zero
+    c.vector()[negs] = 0.0
+    return c
+
+def compute_omt(I0, I1, tau, alpha, space="CG", adjust=False, threshold=1.e-3):
+    """Given two image intensity fields I0 and I1, a time step tau > 0 and
+    a regularization parameter alpha > 0, return a vector field phi
+    such that
+    
+      I_t + div (phi I) \approx 0  in \Omega
+
+    where I(t = 0) = I0 and I(t = tau) = I1. 
+
+    Mueller, M., Karasev, P., Kolesov, I. & Tannenbaum, A. Optical
+    flow estimation for flame detection in videos. IEEE Transactions on
+    image processing, 22, 2786–2797
+
+    If adjust is given as True, filter the input values by replacing
+    negative values by zero, 
+
+    Moreover, set ident_zeros for those rows for which no values
+    exceed the given threshold.
+
+    """
+
+    # Map negative values to zero in I0 and I1:
+    if adjust:
+        info("Adjusting negative input values to zero.")
+        I0 = non_negativize(I0)
+        I1 = non_negativize(I1)
+
+    # Check balance of mass assumption
+    m1 = assemble(I0*dx())
+    m2 = assemble(I1*dx())
+    info("OMT: Mass preserved (\int I0 == \in I1)? I0 = %.5g, I1 = %0.5g" % (m1, m2))
+
+    Q = I0.function_space()
+    mesh = Q.mesh()
+    d = mesh.topology().dim()
+    
+    if space == "RT":
+        V = FunctionSpace(mesh, "RT", 1)
+    else:
+        V = VectorFunctionSpace(mesh, "CG", 1)
+        
+    # Define dI and barI
+    dI = (I1 - I0)/tau
+    barI = 0.5*(I0 + I1)
+    I = barI
+
+    # Define OMT variational forms
+    phi = TrialFunction(V)
+    psi = TestFunction(V)
+    a = (inner(div(I*phi), div(I*psi)) + alpha*inner(I*phi, psi))*dx()
+    L = - inner(dI, div(I*psi))*dx()
+
+    # Assemble linear system
+    A = assemble(a)
+    b = assemble(L)
+
+    # Set zero (tol) rows to identity as relevant
+    A.ident_zeros(tol=threshold)
+
+    # Solve the linear system
+    phi = Function(V)
+    solve(A, phi.vector(), b, "mumps")
+
+    return phi
+
+def element_cell(mesh):
+    "Return the finite element cell."
+    return triangle if mesh.topology().dim() == 2 else tetrahedron
+
+def compute_ocd(c1, c2, tau, D, alpha):
+
+    mesh = c1.function_space().mesh()
+
+    cell = element_cell(mesh)
+    C = FiniteElement("CG", cell, 1) # H1
+    Y = FiniteElement("CG", cell, 1) # H10
+    Q = VectorElement("CG", cell, 1) # H1 (H(div) formulation also natural)
+    Mx = MixedElement([C, Y, Q])
+    M = FunctionSpace(mesh, Mx)
+
+    m = Function(M)
+    (c, y, phi) = split(m)
+    (d, z, psi) = TestFunctions(M)
+
+    F1 = (c*d + (1./tau*d + div(d*phi))*y + inner(D*grad(d), grad(y)) - c2*d)*dx()
+    F2 = (div(c*psi)*y + alpha*(inner(phi, psi) + inner(grad(phi), grad(psi))))*dx()
+    F3 = ((1./tau*c + div(c*phi))*z + inner(D*grad(c), grad(z)) - 1./tau*c1*z)*dx()
+    F = F1 + F2 + F3
+
+    assign(m.sub(0), c2)
+
+    #set_log_level(LogLevel.DEBUG)
+    
+    bcs = [DirichletBC(M.sub(0), c2, "on_boundary"),
+           DirichletBC(M.sub(1), 0.0, "on_boundary")]
+
+    #info(NonlinearVariationalSolver.default_parameters(), True)
+    ps = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "mumps", "maximum_iterations": 50, "relative_tolerance": 1.e-10}}
+    solve(F == 0, m, bcs, solver_parameters=ps)
+    (c, y, phi) = m.split(deepcopy=True)
+
+    return (c, y, phi) 
+
+def compute_magnitude_field(phi):
+
+    # Split phi into its components (really split it)
+    phis = phi.split(deepcopy=True)
+
+    # Get the map from vertex index to dof
+    v2d = [vertex_to_dof_map(phii.function_space()) for phii in phis]
+    dim = len(v2d)
+    
+    # Create array of component values at vertex indices so p0[0] is
+    # value at vertex 0 e.g. 
+    if dim == 2:
+        p0 = phis[0].vector().get_local(v2d[0])
+        p1 = phis[1].vector().get_local(v2d[1])
+        # Take element-wise magnitude
+        m = numpy.sqrt((numpy.square(p0) + numpy.square(p1)))
+    else:
+        p0 = phis[0].vector().get_local(v2d[0])
+        p1 = phis[1].vector().get_local(v2d[1])
+        p2 = phis[2].vector().get_local(v2d[2])
+        # Take element-wise magnitude
+        m = numpy.sqrt((numpy.square(p0) + numpy.square(p1) + numpy.square(p2)))
+        
+    # Map array of values at vertices back to values at dofs and
+    # insert into magnitude field.
+    M = phis[0].function_space()
+    magnitude = Function(M)
+    d2v = dof_to_vertex_map(M)
+    magnitude.vector()[:] = m[d2v]
+    return magnitude
+
+if __name__ == "__main__":
+    mesh = UnitSquareMesh(4, 4)
+    V = VectorFunctionSpace(mesh, "CG", 1)
+    phi = Function(V)
+    phi.vector()[1] = 2.0
+    phi.vector()[0] = 2.0
+    a = compute_magnitude_field(phi)
+    print(a.vector().max())
+    print(a.vector().size())
diff --git a/optimalflow_adjoint.py b/optimalflow_adjoint.py
new file mode 100644
index 0000000..255c6fc
--- /dev/null
+++ b/optimalflow_adjoint.py
@@ -0,0 +1,133 @@
+import resource
+import csv
+import os
+import os.path
+from cases import csvwrite
+
+from dolfin import *
+from dolfin_adjoint import *
+
+def compute_reduced_ocd(c0, c1, tau, D, alpha, results_dir, space="CG", reg="H1"):
+    # c0:
+    # c1:
+    # tau:   time step
+    # D:     diffusion coefficient
+    # alpha: regularization parameter
+    # space: finite element space for the velocity field ("CG" | "RT" | "BDM")
+    # reg:   if "H1", use H1-regularization, else use H(div)
+    
+    info("Computing OCD via reduced approach")
+
+    # Define mesh and function space for the concentration
+    mesh = c0.function_space().mesh()
+    C = FunctionSpace(mesh, "CG", 1)
+    
+    # Space for the convective velocity field phi
+    if space == "CG":
+        Q = VectorFunctionSpace(mesh, "CG", 1)
+    else:
+        Q = FunctionSpace(mesh, space, 1)
+    phi = Function(Q, name="Control")
+
+    # Regularization term
+    def R(phi, alpha, mesh):
+        if reg == "H1":
+            form = 0.5*alpha*(inner(phi, phi) + inner(grad(phi), grad(phi)))*dx(domain=mesh)
+        else:
+            form = 0.5*alpha*(inner(phi, phi) + inner(div(phi), div(phi)))*dx(domain=mesh)
+        return form
+
+    # Define previous solution
+    c_ = Function(C)
+    c_.assign(c0) # Hack to make dolfin-adjoint happy, maybe just start tape here?
+
+    c2 = Function(C)
+    c2.assign(c1)
+    
+    # Define variational problem
+    c = TrialFunction(C)
+    d = TestFunction(C)
+    F = (1.0/tau*(c - c_)*d + div(c*phi)*d + inner(D*grad(c), grad(d)))*dx() 
+    a, L = system(F)
+    bc = DirichletBC(C, c2, "on_boundary")
+
+    # ... and solve it once
+    c = Function(C, name="State")
+    solve(a == L, c, bc, solver_parameters={"linear_solver": "mumps"})
+
+    # Output max values of target and current solution for progress
+    # purposes
+    info("\max c_1 = %f" % c2.vector().max())
+    info("\max c = %f" % c.vector().max())
+
+    # Define the objective functional
+    j = 0.5*(c - c2)**2*dx(domain=mesh) + R(phi, alpha, mesh)
+    J = assemble(j)
+    info("J (initial) = %f" % J)
+    
+    # Define control field
+    m = Control(phi)
+
+    # Define call-back for output at each iteration of the optimization algorithm
+    name = lambda s: os.path.join(results_dir, "opts", s)
+    dirname = os.path.join(results_dir, "opts")
+    if not os.path.isdir(dirname):
+        os.mkdir(dirname)
+    header = ("j", "\max \phi")
+
+    # Make an optimization counter file
+    with open(os.path.join(results_dir, "counter.csv"), "w") as f:
+        writer = csv.writer(f)
+        writer.writerow((0,))     
+        
+    def eval_cb(j, phi):
+        values = (j, phi.vector().max())
+        mem = resource.getrusage(resource.RUSAGE_SELF)[2]
+        info("Current memory usage: %g (MB)" % (mem/1024))
+        info("\tj = %f, \max phi = %f (mm/h)" % values)
+        csvwrite(name("optimization_values.csv"), values, header, mode="a")
+
+        # Read the optimization counter file, update counter, and
+        # write it back, geez. 
+        counter = 0
+        with open(os.path.join(results_dir, "counter.csv"), "r") as f:
+            reader = csv.reader(f)
+            for row in reader:
+                counter = int(row[0])
+        counter += 1
+
+        with open(os.path.join(results_dir, "counter.csv"), "w") as f:
+            info("Updating counter file, counter is now %d " % counter)
+            writer = csv.writer(f)
+            writer.writerow((counter,))     
+
+        # Write current control variable to file in HDF5 and PVD formats
+        file = HDF5File(mesh.mpi_comm(), name("opt_phi_%d.h5" % counter), "w")
+        file.write(phi, "/function", 0)
+        file.close()
+
+    # Define reduced functional in terms of J and m
+    Jhat = ReducedFunctional(J, m, eval_cb_post=eval_cb)
+
+    # Minimize functional
+    tol = 1.0e-12
+    phi_opt = minimize(Jhat,
+                       tol=tol, 
+                       options={"gtol": tol, "maxiter": 500, "disp": True})
+    pause_annotation()
+
+    # Update phi, and do a final solve to compute c
+    phi.assign(phi_opt)
+    solve(a == L, c, bc, solver_parameters={"linear_solver": "mumps"})
+
+    J = assemble(j)
+    j0 = 0.5*(c - c2)**2*dx(domain=mesh)
+    jr = R(phi, alpha, mesh)
+    J0 = assemble(j0)
+    Jr = assemble(jr)
+    info("J  = %f" % J)
+    info("J0 = %f" % J0)
+    info("Jr = %f" % Jr)
+    
+    return (c, phi) 
+
diff --git a/run_case.py b/run_case.py
new file mode 100644
index 0000000..1d482fd
--- /dev/null
+++ b/run_case.py
@@ -0,0 +1,90 @@
+import os
+import os.path
+import pylab
+import numpy
+
+from dolfin import *
+from optimalflow import *
+from cases import *
+
+from optimalflow_adjoint import *
+
+def run_synthetic(n=4, method=None, alpha=1.e-4, D_e=0.0, g=0.0, noise=0.0, output=None,
+                  eps=1.e-4, D=0.0):
+    """ 
+    n:       meshsize 
+    method:  "OMT", "OCD" or "rOCD"
+    alpha:   Regularization parameter (both for OMT, OCD and OCDr)
+    D_e:     Diffusion coefficient for the data generated (and for OCDR)
+    noise:   Noise level for the generated data
+    output:  String with output directory
+
+    eps:     Used for filtering OMT approximation if OMT is given
+    D:       Used as diffusion coefficient for (r)OCD if (r)OCD is given
+    """
+        
+    # Generate synthetic data: c_1, c_2, tau = t_2 - t_1, and exact phi
+    mesh = UnitSquareMesh(n, n)
+    (c1, c2, phi_e, tau) = generate_data_bump(mesh, g, D=D_e, noise=noise, output=output)
+
+    # Compute Linfty, L2 and H10 norms of the exact velocity field phi_e
+    phi_e_mag = compute_magnitude_field(phi_e)
+    phi_e_max = phi_e_mag.vector().max()
+    values = (norm(phi_e, "L2"),  phi_e_max, norm(phi_e, "H10"))
+
+    # Output it
+    filename = lambda f: os.path.join(output, f)
+    name = filename("norms.csv")
+    labels = ("|phi_e|_0", "phi_e_max", "|grad phi_e|_0")
+    csvwrite(name, values, header=labels)
+
+    # Compute the OMT approximation
+    if (method == "OMT"):
+        phi = compute_omt(c1, c2, tau, alpha, space="CG", threshold=eps)
+
+    # or the OCD approximation
+    elif (method == "OCD"):
+        (c, y, phi) = compute_ocd(c1, c2, tau, D, alpha)
+
+    # or the OCD approximation
+    elif (method == "rOCD"):
+        set_working_tape(Tape())
+        (c, phi) = compute_reduced_ocd(c1, c2, tau, D, alpha, output)
+
+    # Store computed solution to .h5 and .pvd format
+    hdf5write(filename("phi.h5"), phi)
+    file = File(filename("phi.pvd"))
+    file << phi
+
+    # Also store the concentration c if we estimate it
+    if not (method == "OMT"):
+        hdf5write(filename("c.h5"), c)
+        file = File(filename("c.pvd"))
+        file << c
+    
+    # Compute magnitude of the field, and extract peak magnitude
+    m = compute_magnitude_field(phi)
+    phi_max = m.vector().max()
+        
+    labels = ("|phi|_0", "phi_max", "|grad phi|_0")
+    values = (norm(phi, "L2"),  phi_max, norm(phi, "H10"))
+    csvwrite(name, values, header=labels, mode="a")
+
+if __name__ == "__main__":
+
+    # Turn on FEniCS optimizations
+    parameters["form_compiler"]["cpp_optimize"] = True
+    flags = ["-O3", "-ffast-math", "-march=native"]
+    parameters["form_compiler"]["cpp_optimize_flags"] = " ".join(flags)
+
+    n = 32
+    g = 2.e-2                # Velocity field phi = (g, g)
+    D = 1.e-2                # Diffusion coefficient (for data and OCD)
+    noise = 0                # Data generation noise level
+    method= "rOCD"           # "OMT", "OCD" or "rOCD"
+    alpha = 1.e-3            # Regularization parameter
+    eps = 1.e-6              # OMT threshold
+    results = os.path.join("results-r0", "method_%s_n_%d_alpha_%2.1g" % (method, n, alpha))
+    
+    run_synthetic(n=n, method=method, alpha=alpha, D_e=D, g=g, noise=noise,
+                  output=results, eps=eps, D=D)