diff --git a/scripts/run_experiments.jl b/scripts/run_experiments.jl
index 170bed8166f77a0eee807aff131abbe1e06aed45..954c1935c72e48cc823ee4e687d802d8bf6ec222 100644
--- a/scripts/run_experiments.jl
+++ b/scripts/run_experiments.jl
@@ -794,10 +794,14 @@ function denoise_pd(ctx)
 
     #project_l2_lagrange!(st.g, ctx.params.g_arr)
     interpolate!(st.g, x -> evaluate_bilinear(ctx.params.g_arr, x))
-    st.u.data .= st.g.data
+    if haskey(ctx.params, :init_zero) && ctx.params.init_zero
+        st.u.data .= 0
+    else
+        st.u.data .= st.g.data
+    end
 
-    save_denoise(st, i) =
-	output(st, "output/$(st.name)_$(lpad(i, 5, '0')).vtu",
+    save_step(st, i) =
+        output(st, joinpath(ctx.outdir, "output_$(ctx.params.algorithm)_$(lpad(i, 5, '0')).vtu"),
 	    st.g, st.u, st.p1, st.p2)
 
     log!(x::Nothing; kwargs...) = x
@@ -808,42 +812,41 @@ function denoise_pd(ctx)
 	println(row)
     end
 
-    #pvd = paraview_collection("output/$(st.name).pvd")
-    #pvd[0] = save_denoise(st, 0)
+    pvd = paraview_collection(joinpath(ctx.outdir, "output_$(ctx.params.algorithm).pvd")) do pvd
+        k = 0
+        println("primal energy: $(primal_energy(st))")
 
-    k = 0
-    println("primal energy: $(primal_energy(st))")
+        while true
+            k += 1
+            if ctx.params.algorithm == :pd1
+                # no step size control
+                step_pd2!(st; sigma, tau, theta)
+            elseif ctx.params.algorithm == :pd2
+                theta = 1 / sqrt(1 + 2 * mu * tau)
+                tau *= theta
+                sigma /= theta
+                step_pd2!(st; sigma, tau, theta)
+            elseif ctx.params.algorithm == :newton
+                step!(st)
+            end
+            #pvd[k] = save_step(st, k)
 
-    while true
-	k += 1
-	if ctx.params.algorithm == :pd1
-	    # no step size control
-	    step_pd2!(st; sigma, tau, theta)
-	elseif ctx.params.algorithm == :pd2
-	    theta = 1 / sqrt(1 + 2 * mu * tau)
-	    tau *= theta
-	    sigma /= theta
-	    step_pd2!(st; sigma, tau, theta)
-	elseif ctx.params.algorithm == :newton
-	    step!(st)
-	end
-	#pvd[k] = save_denoise(st, k)
+            domain_factor = 1 / sqrt(area(mesh))
+            norm_step_ = norm_step(st) * domain_factor
+            #estimate_pd!(st)
 
-	domain_factor = 1 / sqrt(area(mesh))
-	norm_step_ = norm_step(st) * domain_factor
-        #estimate_pd!(st)
+            log!(ctx.params.df; k,
+                norm_step = norm_step_,
+                #est = estimate_error(st),
+                primal_energy = primal_energy(st))
 
-        log!(ctx.params.df; k,
-            norm_step = norm_step_,
-            #est = estimate_error(st),
-            primal_energy = primal_energy(st))
+            #norm_residual_ < params.tol && norm_step_ < params.tol && break
+            haskey(params, :tol) && norm_step_ < params.tol && break
+            haskey(params, :max_iters) && k >= params.max_iters && break
+        end
 
-	#norm_residual_ < params.tol && norm_step_ < params.tol && break
-	haskey(params, :tol) && norm_step_ < params.tol && break
-        haskey(params, :max_iters) && k >= params.max_iters && break
+        pvd[1] = save_step(st, 0)
     end
-    #pvd[1] = save_denoise(st, 1)
-    #vtk_save(pvd)
     return st
 end
 
@@ -896,6 +899,66 @@ function experiment_convergence_rate(ctx)
         energy_min, algparams...)
 end
 
+function experiment_convergence_difference(ctx)
+    g_arr = from_img([0. 0.; 1. 0.])
+    mesh = init_grid(g_arr;)
+
+    algparams = (
+        init_zero = true,
+        alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
+        gamma1 = 1e-10, gamma2 = 1e-10,
+        tol = 1e-10,
+        max_iters = 10000,
+    )
+    algctx = Util.Context(ctx; g_arr, mesh, algparams...)
+
+    df1 = DataFrame()
+    df2 = DataFrame()
+    st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
+    st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
+    p1 = plot(df1.k, df1.norm_step; axis = :log)
+    p2 = plot(df2.k, df2.norm_step; axis = :log)
+    p3 = plot(df1.k, df1.primal_energy; xaxis = :log)
+    p4 = plot(df2.k, df2.primal_energy; xaxis = :log)
+    display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
+
+    display(st1.u.data)
+    display(st2.u.data)
+
+    return
+end
+
+function experiment_convergence_oscillations(ctx)
+    g_arr = from_img([0. 0.; 1. 0.])
+    mesh = init_grid(g_arr;)
+
+    algparams = (
+        #alpha1 = 0., alpha2 = 30., lambda = 1., beta = 0.,
+        #alpha1 = 0.5, alpha2 = 0., lambda = 0., beta = 1.,
+        alpha1 = 1., alpha2 = 0., lambda = 0., beta = 1., # oscillating
+        # DP0 for p1 has less oscillations
+        gamma1 = 1e-5, gamma2 = 1e-5,
+        tol = 1e-10,
+        max_iters = 10000,
+    )
+    algctx = Util.Context(ctx; g_arr, mesh, algparams...)
+
+    df1 = DataFrame()
+    df2 = DataFrame()
+    st1 = denoise_pd(Util.Context(algctx; algorithm = :pd2, df = df1));
+    st2 = denoise_pd(Util.Context(algctx; algorithm = :newton, df = df2));
+    p1 = plot(df1.k, df1.norm_step; axis = :log)
+    p2 = plot(df2.k, df2.norm_step; axis = :log)
+    p3 = plot(df1.k, df1.primal_energy; xaxis = :log)
+    p4 = plot(df2.k, df2.primal_energy; xaxis = :log)
+    display(plot(p1, p2, p3, p4; layout = @layout[a b; c d]))
+
+    display(st1.u.data)
+    display(st2.u.data)
+
+    return
+end
+
 function denoise_approximation(ctx)
     domain_factor = 1 / sqrt(area(ctx.params.mesh))
 
@@ -1016,6 +1079,7 @@ function denoise_approximation(ctx)
     return st
 end
 
+
 function experiment_approximation(ctx)
     #g_arr = from_img([1. 0.; 0. 0.])
     g_arr = from_img([0. 0.; 1. 0.])