diff --git a/.gitignore b/.gitignore
index 8d69778b36a7e14165c6bceba7b6a0bc6a0edeb9..1fcac93c31da5514aeaed40deef9a86ce4713775 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
 *.vtu
+*.pvd
diff --git a/src/function.jl b/src/function.jl
index 9cbd6e65b5b49b3803d5a5a9b3c90de151c79e1d..e3e3cf0dbbb20572b7622ea86f33c8cc8ccede73 100644
--- a/src/function.jl
+++ b/src/function.jl
@@ -183,19 +183,22 @@ end
 
 
 
-append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.space.element)
+vtk_append!(vtkfile, f::FeFunction, fs::FeFunction...) =
+    (vtk_append!(vtkfile, f); vtk_append!(vtkfile, fs...))
 
-function append_data!(vtkfile, f::FeFunction, ::P1)
-    # separate out data per cell
+vtk_append!(vtkfile, f::FeFunction) =
+    vtk_append!(vtkfile, f, f.space.element)
+
+function vtk_append!(vtkfile, f::FeFunction, ::P1)
+    # separate out data per cell since different cells have distinct vertices
+    # in the output vtk
     fd = vec(f.data[f.space.dofmap])
     vtk_point_data(vtkfile, fd, f.name)
 end
 
-function append_data!(vtkfile, f::FeFunction, ::DP0)
+function vtk_append!(vtkfile, f::FeFunction, ::DP0)
     vtk_cell_data(vtkfile, f.data, f.name)
 end
 
-append_data!(vtkfile, f::FeFunction, ::DP1) =
-    append_data!(vtkfile, f, P1())
-
-
+vtk_append!(vtkfile, f::FeFunction, ::DP1) =
+    vtk_append!(vtkfile, f, P1())
diff --git a/src/mesh.jl b/src/mesh.jl
index 3835b17d1b4aa1fd4ea62cf395e6a2d4cffac8a3..02b2d8f7f5c80507fa5ecf8991fa9465be9d550d 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -43,17 +43,21 @@ init_grid(img::Array{<:Any, 2}, type=:vertex) =
 #    return all(λ .>= 0) && sum(λ) <= 1
 #end
 
-save(filename, f::FeFunction, fs::FeFunction...) =
-    save(filename, f.space.mesh, f, fs...)
+function vtk_mesh(filename, mesh::Mesh)
+    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, 3*(i-1)+1:3*(i-1)+3)
+	for i in axes(mesh.cells, 2)]
+    vertices = reshape(mesh.vertices[:, mesh.cells], 2, :)
+
+    return vtk_grid(filename, vertices, cells)
+end
 
+"convenience function for saving to vtk"
 function save(filename, mesh::Mesh, fs...)
-    cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, 3*(i-1)+1:3*(i-1)+3) for i in axes(mesh.cells, 2)]
-    vertices = reshape(mesh.vertices[:, mesh.cells], 2, :)
-    vtkfile = vtk_grid(filename, vertices, cells)
+    vtk = vtk_mesh(filename, mesh)
     for f in fs
 	f.space.mesh == mesh ||
 	    throw(ArgumentError("meshes do not match"))
-        append_data!(vtkfile, f)
+        vtk_append!(vtk, f)
     end
-    vtk_save(vtkfile)
+    vtk_save(vtk)
 end
diff --git a/src/run.jl b/src/run.jl
index 47a49c7b4717e2a6529d546e197d766dcf59c70d..6b81628faf4a9b8d9f8b30a5865e34de3768c368 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -5,6 +5,7 @@ using LinearAlgebra: norm
 
 
 function myrun()
+    name = "test"
     mesh = init_grid(50)
     d = 2
     m = 1
@@ -63,6 +64,21 @@ function myrun()
 	return -p2 + lambda / m2 * (nablau + nabladu) - cond
     end
 
+    function p1_project!()
+	p1.space.element::DP0
+	p1.data .= clamp.(p1.data, -alpha1, alpha1)
+    end
+
+    function p2_project!()
+	p2.space.element::DP1
+	p2d = reshape(p2.data, m * d * 3, :) # no copy
+	for i in axes(p2d, 2)
+	    p2in = norm(p2d[:, i])
+	    if p2in > lambda
+		p2d[:, i] .*= lambda ./ p2in
+	    end
+	end
+    end
 
     @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2)
 	m1 = max(gamma1, norm(T(u) - g))
@@ -97,31 +113,35 @@ function myrun()
 	return -aB - p1part - p2part + gpart
     end
 
+    function save_step!(pvd, i)
+	vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh)
+	vtk_append!(vtk, g, u, p1, p2)
+	vtk_save(vtk)
+	pvd[i] = vtk
+    end
+
+    pvd = paraview_collection("$(name).pvd")
+    save_step!(pvd, 0)
     for i = 1:6
+	# solve du
 	A = assemble(Vu, du_a; g, u, nablau, p1, p2)
 	b = assemble_rhs(Vu, du_l; g, u, nablau)
-
 	du.data .= A \ b
 
+	# solve dp1, dp2
 	interpolate!(dp1, dp1_update; g, u, p1, du)
 	interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
 
+	# newton update
 	u.data .+= du.data
 	p1.data .+= dp1.data
 	p2.data .+= dp2.data
 
-	# reproject
-	p1.space.element::DP0
-	p1.data .= clamp.(p1.data, -alpha1, alpha1)
-	p2.space.element::DP1
-	p2d = reshape(p2.data, m * d * 3, :) # no copy
-	for i in axes(p2d, 2)
-	    p2in = norm(p2d[:, i])
-	    if p2in > lambda
-		p2d[:, i] .*= lambda ./ p2in
-	    end
-	end
-    end
+	# reproject p1, p2
+	p1_project!()
+	p2_project!()
 
-    save("test.vtu", g, u, p1, p2)
+	save_step!(pvd, i)
+    end
+    vtk_save(pvd)
 end