diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index bc3f44708ad77b2353933901cec74cf99ff378d8..321fee644ea5b6b1a0d7aab0f2567af6b256799f 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -77,7 +77,7 @@ function L1L2TVContext(name, mesh, m; T, tdata, S,
     Vest = FeSpace(mesh, DP0(), (1,))
     Vg = FeSpace(mesh, P1(), (1,))
     Vu = FeSpace(mesh, P1(), (m,))
-    Vp1 = FeSpace(mesh, DP1(), (1,))
+    Vp1 = FeSpace(mesh, P1(), (1,))
     Vp2 = FeSpace(mesh, DP0(), (m, d))
 
     est = FeFunction(Vest, name="est")
@@ -399,6 +399,7 @@ function estimate!(ctx::L1L2TVContext)
 
     w = FeFunction(ctx.u.space)
     solve_primal!(w, ctx)
+    #w.data .= .-w.data
     # TODO: find better name: is actually a cell-wise integration
     project!(ctx.est, estf; ctx.g, ctx.u, ctx.p1, ctx.p2,
         nablau = nabla(ctx.u), w, nablaw = nabla(w), ctx.tdata)
@@ -679,6 +680,7 @@ function denoise_approximation(ctx)
 
     T(tdata, u) = u
     S(u, nablau) = u
+    #S(u, nablau) = nablau
 
     st = L1L2TVContext(ctx.params.name, ctx.params.mesh, 1;
         T, tdata = nothing, S,
@@ -688,6 +690,7 @@ function denoise_approximation(ctx)
     #project_img!(st.g, img)
     interpolate!(st.g, x -> interpolate_bilinear(ctx.params.img, x))
     st.u.data .= st.g.data
+    #st.u.data .= rand(size(st.u.data))
 
     save_vtk(st, i) =
         output(st,
@@ -708,6 +711,24 @@ function denoise_approximation(ctx)
         pvd[0] = save_vtk(st, 0)
         println("primal energy: $(primal_energy(st))")
 
+        # initial refinement
+        for _ in 1:14
+            marked_cells = Set(axes(st.mesh.cells, 2))
+            mesh, fs = refine(st.mesh, marked_cells;
+                st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
+            st = L1L2TVContext("run", mesh, st.d, st.m, T, nothing, 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)
+        end
+        #marked_cells = Set(axes(st.mesh.cells, 2))
+        #mesh2, fs = refine(st.mesh, marked_cells;
+        #    st.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2)
+        #st2 = L1L2TVContext("run", mesh2, st.d, st.m, T, nothing, 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)
+
         k = 0
         while true
             step!(st)