diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 93ccb9d533715e907da4a067f06cafaba4589a5c..5bbf6eb01a815b30416e65409ff0723a5964b7e0 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -851,18 +851,19 @@ function inpaint(img, imgmask; name, params...)
 end
 
 function optflow(ctx)
-    tol_newton = 1e-6 # cauchy criterion
-
     size(ctx.params.imgf0) == size(ctx.params.imgf1) ||
 	throw(ArgumentError("non-matching image sizes"))
 
+    tol_newton = 1e-5 # cauchy criterion for inner newton loop
+
     m = 2 # range dim of vector field
+
+    # convert to cartesian coordinates
     imgf0 = from_img(ctx.params.imgf0)
     imgf1 = from_img(ctx.params.imgf1)
 
-    #mesh = init_grid(imgf0)
-    mesh = init_grid(imgf0, (size(imgf0) .÷ 2^3)...)
-    domain_factor = 1 / sqrt(area(mesh))
+    mesh = init_grid(imgf0, (size(imgf0) .÷ 2^2)...)
+    mesh_area = area(mesh)
 
     # optflow specific stuff
     Vg = FeSpace(mesh, P1(), (1,))
@@ -872,19 +873,21 @@ function optflow(ctx)
 
     T(tdata, u) = tdata * u # here tdata will be nabla(fw)
     S(u, nablau) = nablau
-    #S(u, nablau) = u
 
     st = L1L2TVContext("run", mesh, m; T, tdata = nabla(fw), S,
         ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
+    u_acc = FeFunction(st.u.space, name="u_acc")
+    u_acc.data .= 0
+
     function warp!()
-        imgfw = warp_backwards(imgf1, sample(st.u))
-        project_l2_lagrange!(fw, imgfw)
+        imgfw = warp_backwards(imgf1, sample(u_acc))
+        project_l2_pixel!(fw, imgfw)
 
         g_optflow(x; u, f0, fw, nablafw) =
-            #nablafw * u
-            nablafw * u - (fw - f0)
+            -(fw - f0)
+            #nablafw * u - (fw - f0)
         interpolate!(st.g, g_optflow; st.u, f0, fw, nablafw = st.tdata)
     end
 
@@ -898,7 +901,7 @@ function optflow(ctx)
 
     save_step(i) =
         output(st, joinpath(ctx.outdir, "output_$(lpad(i, 5, '0')).vtu"),
-	    st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw)
+	    st.g, st.u, st.p1, st.p2, st.est, f0, f1, fw, u_acc)
 
     i = 0
     pvd = paraview_collection(joinpath(ctx.outdir, "output.pvd")) do pvd
@@ -917,34 +920,38 @@ function optflow(ctx)
                 #estimate!(st)
                 pvd[i] = save_step(i)
 
-                norm_step_ = norm_step(st) * domain_factor
+                norm_step_ = norm_step(st) / sqrt(mesh_area)
                 println("norm_step = $norm_step_")
 
                 # residual is basically unusable: takes long to compute and
                 # oscillates
-                #norm_residual_ = norm_residual(st) * domain_factor
+                #norm_residual_ = norm_residual(st) / sqrt(mesh_area)
                 #println("norm_residual = $norm_residual_")
 
-                display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
+                #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
             end
+
+            u_acc.data .+= st.u.data
+            display(plot(colorflow(to_img(sample(u_acc)); ctx.params.maxflow)))
+
             # poor man's hash
             #println(integrate(st.mesh, (x; tdata) -> tdata; st.tdata))
 
-            println("refine ...")
-            marked_cells = Set(axes(mesh.cells, 2))
-            mesh, fs = refine(mesh, marked_cells;
-                st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2,
-                f0, f1, fw)
-            st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fs.fw), S,
-                st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2,
-                fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
-            f0, f1, fw = (fs.f0, fs.f1, fs.fw)
+            #println("refine ...")
+            #marked_cells = Set(axes(mesh.cells, 2))
+            #mesh, fs = refine(mesh, marked_cells;
+            #    st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2,
+            #    f0, f1, fw, u_acc)
+            #st = L1L2TVContext("run", mesh, st.d, st.m, T, nabla(fs.fw), S,
+            #    st.alpha1, st.alpha2, st.beta, st.lambda, st.gamma1, st.gamma2,
+            #    fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
+            #f0, f1, fw, u_acc = (fs.f0, fs.f1, fs.fw, fs.u_acc)
 
-            #i += 1
-            pvd[i] = save_step(i)
+            ##i += 1
+            #pvd[i] = save_step(i)
 
-            println("interpolate image data ...")
-            interpolate_image_data!()
+            #println("interpolate image data ...")
+            #interpolate_image_data!()
 
             println("warp ...")
             warp!() # yay
@@ -978,7 +985,7 @@ function optflow(ctx)
         #    i += 1
         #    pvd[i] = save_step(i)
     end
-    display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
+    #display(plot(colorflow(to_img(sample(st.u)); ctx.params.maxflow)))
 
     #CSV.write(joinpath(ctx.outdir, "energies.csv"), df)
     #saveimg(joinpath(ctx.outdir, "output_glob.png"), fetch_u(states.glob))