diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 1a9e06587f18026a90283ad5881757c6382b7318..d288ef042da62b86f8a8cc54f9a0976f6fc3ce1b 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -764,15 +764,19 @@ function denoise_approximation(ctx)
         ctx.params.alpha1, ctx.params.alpha2, ctx.params.lambda, ctx.params.beta,
         ctx.params.gamma1, ctx.params.gamma2)
 
+    # primal reconstruction
+    w = FeFunction(st.u.space, name = "w")
+    w.data .= 0
+
     #project_l2_lagrange!(st.g, img)
-    interpolate!(st.g, x -> evaluate_bilinear(ctx.params.img, x))
+    interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
     st.u.data .= st.g.data
     #st.u.data .= rand(size(st.u.data))
 
     save_vtk(st, i) =
         output(st,
             joinpath(ctx.outdir, "$(ctx.params.name)_$(lpad(i, 5, '0')).vtu"),
-	    st.g, st.u, st.p1, st.p2, st.est)
+	    st.g, st.u, st.p1, st.p2, st.est, w)
 
     log!(x::Nothing; kwargs...) = x
     function log!(df::DataFrame; kwargs...)
@@ -789,14 +793,15 @@ function denoise_approximation(ctx)
         println("primal energy: $(primal_energy(st))")
 
         # initial refinement
-        for _ in 1:14
+        for _ in 1:0
             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.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2, w)
             st = L1L2TVState(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)
+            w = fs.w
         end
         #marked_cells = Set(axes(st.mesh.cells, 2))
         #mesh2, fs = refine(st.mesh, marked_cells;
@@ -807,8 +812,9 @@ function denoise_approximation(ctx)
         #    fs.est, fs.g, fs.u, fs.p1, fs.p2, fs.du, fs.dp1, fs.dp2)
 
         k = 0
-        while true
+        while true && k < 50
             step!(st)
+            solve_primal!(w, st)
 
             norm_step_ = norm_step(st) * domain_factor
 
@@ -844,11 +850,12 @@ function denoise_approximation(ctx)
                     mark(st; theta = 0.5) : # estimator + dörfler
                     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.est, st.g, st.u, st.p1, st.p2, st.du, st.dp1, st.dp2, w)
                 st = L1L2TVState(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)
+                w = fs.w
 
                 norm_residual_ = norm_residual(st) * domain_factor
                 estimate!(st)
@@ -870,18 +877,18 @@ function denoise_approximation(ctx)
 end
 
 function experiment_approximation(ctx)
-    df = DataFrame()
-    #img = from_img([1. 0.; 0. 0.])
-    img = from_img([0. 0.; 1. 0.])
-    #mesh = init_grid(img; type=:vertex)
-    mesh = init_grid(img;)
+    #g_arr = from_img([1. 0.; 0. 0.])
+    g_arr = from_img([0. 0.; 1. 0.])
+    mesh = init_grid(g_arr;)
 
+    df = DataFrame()
     denoise_approximation(Util.Context(ctx; name = "test", df,
-        img, mesh,
-        alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
+        g_arr, mesh,
+        #alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
         #alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
-        gamma1 = 1e-4, gamma2 = 1e-4,
-        tol = 1e-10, adaptive = false,
+        alpha1 = 1., alpha2 = 0., lambda = 1., beta = 1e-5,
+        gamma1 = 1e-3, gamma2 = 1e-3,
+        tol = 1e-10, adaptive = true,
     ))
 end