diff --git a/docs/src/reference/io.md b/docs/src/reference/io.md
index cc6fccc6..5adc2248 100644
--- a/docs/src/reference/io.md
+++ b/docs/src/reference/io.md
@@ -1,14 +1,10 @@
 # IO
 
-## Optimiser history and visualisation
+## Optimiser history
 ```@docs
 LevelSetTopOpt.write_history
 ```
 
-```@docs
-LevelSetTopOpt.write_vtk
-```
-
 ## Object IO in serial
 ```@docs
 LevelSetTopOpt.save
diff --git a/scripts/MPI/3d_elastic_compliance_ALM.jl b/scripts/MPI/3d_elastic_compliance_ALM.jl
index 3e5f6f24..8ae4b8c3 100644
--- a/scripts/MPI/3d_elastic_compliance_ALM.jl
+++ b/scripts/MPI/3d_elastic_compliance_ALM.jl
@@ -28,7 +28,9 @@ function main(mesh_partition,distribute,el_size)
   η_coeff = 2
   α_coeff = 4
   vf = 0.4
-  path = dirname(dirname(@__DIR__))*"/results/3d_elastic_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_elastic_compliance_ALM/"
+  iter_mod = 10
+  i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
   model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size);
@@ -98,15 +100,15 @@ function main(mesh_partition,distribute,el_size)
   )
 
   ## Optimiser
-  i_am_main(ranks) && mkdir(path)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_hyperelastic_compliance_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_ALM.jl
index 7bb9a0b2..ad197331 100644
--- a/scripts/MPI/3d_hyperelastic_compliance_ALM.jl
+++ b/scripts/MPI/3d_hyperelastic_compliance_ALM.jl
@@ -24,7 +24,9 @@ function main(mesh_partition,distribute,el_size)
   η_coeff = 2
   α_coeff = 4
   vf = 0.4
-  path = dirname(dirname(@__DIR__))*"/results/3d_hyperelastic_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_hyperelastic_compliance_ALM/"
+  iter_mod = 10
+  i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
   model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size)
@@ -104,15 +106,15 @@ function main(mesh_partition,distribute,el_size)
   )
   
   ## Optimiser
-  i_am_main(ranks) && mkdir(path)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl
index c394a268..0ddd1da4 100644
--- a/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl
+++ b/scripts/MPI/3d_hyperelastic_compliance_neohook_ALM.jl
@@ -24,7 +24,9 @@ function main(mesh_partition,distribute,el_size)
   η_coeff = 2
   α_coeff = 4
   vf = 0.4
-  path = dirname(dirname(@__DIR__))*"/results/3d_hyperelastic_compliance_neohook_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_hyperelastic_compliance_neohook_ALM/"
+  iter_mod = 10
+  i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
   model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size)
@@ -118,15 +120,15 @@ function main(mesh_partition,distribute,el_size)
   )
   
   ## Optimiser
-  i_am_main(ranks) && mkdir(path)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_inverse_homenisation_ALM.jl b/scripts/MPI/3d_inverse_homenisation_ALM.jl
index 488ca57e..2769375b 100644
--- a/scripts/MPI/3d_inverse_homenisation_ALM.jl
+++ b/scripts/MPI/3d_inverse_homenisation_ALM.jl
@@ -26,8 +26,9 @@ function main(mesh_partition,distribute,el_size)
   η_coeff = 2
   α_coeff = 4
   vf = 0.5
-  path = dirname(dirname(@__DIR__))*"/results/3d_inverse_homenisation_ALM"
-  i_am_main(ranks) && mkdir(path)
+  path = dirname(dirname(@__DIR__))*"/results/3d_inverse_homenisation_ALM/"
+  iter_mod = 10
+  i_am_main(ranks) && mkdir(path)  
 
   ## FE Setup
   model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size,isperiodic=(true,true,true))
@@ -106,11 +107,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_inverter_ALM.jl b/scripts/MPI/3d_inverter_ALM.jl
index 4fff7765..1e34d67e 100644
--- a/scripts/MPI/3d_inverter_ALM.jl
+++ b/scripts/MPI/3d_inverter_ALM.jl
@@ -32,7 +32,8 @@ function main(mesh_partition,distribute,el_size)
   δₓ=0.5
   ks = 0.01
   g = VectorValue(1,0,0)
-  path = dirname(dirname(@__DIR__))*"/results/3d_inverter_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_inverter_ALM/"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -120,11 +121,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_inverter_HPM.jl b/scripts/MPI/3d_inverter_HPM.jl
index ad041a18..1e051812 100644
--- a/scripts/MPI/3d_inverter_HPM.jl
+++ b/scripts/MPI/3d_inverter_HPM.jl
@@ -32,7 +32,8 @@ function main(mesh_partition,distribute,el_size)
   δₓ=0.5
   ks = 0.01
   g = VectorValue(1,0,0)
-  path = dirname(dirname(@__DIR__))*"/results/3d_inverter_HPM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_inverter_HPM/"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -121,11 +122,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;γ,γ_reinit,ls_enabled,α_min=0.7,
     ls_γ_max=0.05,verbose=i_am_main(ranks),constraint_names=[:Vol,:UΓ_out])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl b/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl
index 2fb8ce00..19598021 100644
--- a/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl
+++ b/scripts/MPI/3d_nonlinear_thermal_compliance_ALM.jl
@@ -31,7 +31,8 @@ function main(mesh_partition,distribute,el_size)
   η_coeff = 2
   α_coeff = 4
   vf = 0.4
-  path = dirname(dirname(@__DIR__))*"/results/3d_nonlinear_thermal_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/3d_nonlinear_thermal_compliance_ALM/"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -108,11 +109,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl
index cae68f2d..d9da92ff 100644
--- a/scripts/MPI/3d_thermal_compliance_ALM.jl
+++ b/scripts/MPI/3d_thermal_compliance_ALM.jl
@@ -38,6 +38,7 @@ function main(mesh_partition,distribute,el_size,coef,step)
   α_coeff = 4
   vf = 0.4
   path = dirname(dirname(@__DIR__))*"/results/3d_thermal_compliance_ALM_Nx$(el_size[1])_Coef$(coef)_Step$(step)"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -111,11 +112,12 @@ function main(mesh_partition,distribute,el_size,coef,step)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/elastic_compliance_ALM.jl b/scripts/MPI/elastic_compliance_ALM.jl
index 5be16714..49f97724 100644
--- a/scripts/MPI/elastic_compliance_ALM.jl
+++ b/scripts/MPI/elastic_compliance_ALM.jl
@@ -29,6 +29,7 @@ function main(mesh_partition,distribute,el_size)
   α_coeff = 4
   vf = 0.4
   path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_ALM_MPI"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -100,11 +101,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl
index e5399ae1..4b866b5f 100644
--- a/scripts/MPI/inverse_homenisation_ALM.jl
+++ b/scripts/MPI/inverse_homenisation_ALM.jl
@@ -27,6 +27,7 @@ function main(mesh_partition,distribute,el_size)
   α_coeff = 4
   vf = 0.5
   path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM_MPI"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -101,11 +102,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/thermal_compliance_ALM.jl b/scripts/MPI/thermal_compliance_ALM.jl
index f4208c11..b4f316e2 100644
--- a/scripts/MPI/thermal_compliance_ALM.jl
+++ b/scripts/MPI/thermal_compliance_ALM.jl
@@ -29,6 +29,7 @@ function main(mesh_partition,distribute)
   α_coeff = 4
   vf = 0.4
   path = dirname(dirname(@__DIR__))*"/results/thermal_compliance_ALM_MPI"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -86,11 +87,12 @@ function main(mesh_partition,distribute)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,
     verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/MPI/thermal_compliance_ALM_PETSc.jl b/scripts/MPI/thermal_compliance_ALM_PETSc.jl
index a39f1a6a..8c032c04 100644
--- a/scripts/MPI/thermal_compliance_ALM_PETSc.jl
+++ b/scripts/MPI/thermal_compliance_ALM_PETSc.jl
@@ -29,6 +29,7 @@ function main(mesh_partition,distribute,el_size)
   α_coeff = 4
   vf = 0.4
   path = dirname(dirname(@__DIR__))*"/results/thermal_compliance_ALM_MPI+PETSc"
+  iter_mod = 10
   i_am_main(ranks) && mkdir(path)
 
   ## FE Setup
@@ -100,11 +101,12 @@ function main(mesh_partition,distribute,el_size)
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,
     verbose=i_am_main(ranks),constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history;ranks=ranks)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/Serial/elastic_compliance_ALM.jl b/scripts/Serial/elastic_compliance_ALM.jl
index a193d39b..c1f780db 100644
--- a/scripts/Serial/elastic_compliance_ALM.jl
+++ b/scripts/Serial/elastic_compliance_ALM.jl
@@ -26,7 +26,8 @@ function main()
   α_coeff = 4
   vf = 0.4
   g = VectorValue(0,-1)
-  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -85,11 +86,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/elastic_compliance_HPM.jl b/scripts/Serial/elastic_compliance_HPM.jl
index 825e92f9..4e620db2 100644
--- a/scripts/Serial/elastic_compliance_HPM.jl
+++ b/scripts/Serial/elastic_compliance_HPM.jl
@@ -26,7 +26,8 @@ function main()
   α_coeff = 4
   vf = 0.4
   g = VectorValue(0,-1)
-  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_HPM"
+  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_HPM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -86,11 +87,12 @@ function main()
   optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main();
\ No newline at end of file
diff --git a/scripts/Serial/elastic_compliance_LM.jl b/scripts/Serial/elastic_compliance_LM.jl
index 64b04677..b1be54a6 100644
--- a/scripts/Serial/elastic_compliance_LM.jl
+++ b/scripts/Serial/elastic_compliance_LM.jl
@@ -24,7 +24,8 @@ function main()
   η_coeff = 2
   α_coeff = 4
   g = VectorValue(0,-1)
-  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_LM"
+  path = dirname(dirname(@__DIR__))*"/results/elastic_compliance_LM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -81,11 +82,12 @@ function main()
   ## Optimiser
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;γ,γ_reinit,verbose=true)
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main();
\ No newline at end of file
diff --git a/scripts/Serial/hyperelastic_compliance_ALM.jl b/scripts/Serial/hyperelastic_compliance_ALM.jl
index f2117446..e6b801dc 100644
--- a/scripts/Serial/hyperelastic_compliance_ALM.jl
+++ b/scripts/Serial/hyperelastic_compliance_ALM.jl
@@ -17,7 +17,8 @@ function main()
   η_coeff = 2
   α_coeff = 4
   vf = 0.6
-  path = dirname(dirname(@__DIR__))*"/results/hyperelastic_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/hyperelastic_compliance_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -85,11 +86,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
index 0ebae329..1a10577c 100644
--- a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
+++ b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
@@ -17,7 +17,8 @@ function main()
   η_coeff = 2
   α_coeff = 4
   vf = 0.6
-  path = dirname(dirname(@__DIR__))*"/results/hyperelastic_compliance_neohook_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/hyperelastic_compliance_neohook_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -99,11 +100,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl
index d4d60aa8..7ac70885 100644
--- a/scripts/Serial/inverse_homenisation_ALM.jl
+++ b/scripts/Serial/inverse_homenisation_ALM.jl
@@ -24,7 +24,8 @@ function main()
   η_coeff = 2
   α_coeff = 4
   vf = 0.5
-  path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/inverse_homenisation_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -85,11 +86,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/inverter_ALM.jl b/scripts/Serial/inverter_ALM.jl
index 85b28d89..9b1e4879 100644
--- a/scripts/Serial/inverter_ALM.jl
+++ b/scripts/Serial/inverter_ALM.jl
@@ -30,7 +30,8 @@ function main()
   δₓ = 0.2
   ks = 0.1
   g = VectorValue(0.5,0)
-  path = dirname(dirname(@__DIR__))*"/results/inverter_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/inverter_ALM/"
+  iter_mod = 10
   mkdir(path)
   
   ## FE Setup
@@ -100,11 +101,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol,:UΓ_out])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/inverter_HPM.jl b/scripts/Serial/inverter_HPM.jl
index f79ca7e9..4f22c0f3 100644
--- a/scripts/Serial/inverter_HPM.jl
+++ b/scripts/Serial/inverter_HPM.jl
@@ -30,7 +30,8 @@ function main()
   δₓ = 0.2
   ks = 0.1
   g = VectorValue(0.5,0)
-  path = dirname(dirname(@__DIR__))*"/results/inverter_HPM_osc"
+  path = dirname(dirname(@__DIR__))*"/results/inverter_HPM/"
+  iter_mod = 10
   mkdir(path)
   
   ## FE Setup
@@ -102,11 +103,12 @@ function main()
   optimiser = HilbertianProjection(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,α_min=0.4,ls_enabled,verbose=true,constraint_names=[:Vol,:UΓ_out])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
index 270eedcf..ba3f8d74 100644
--- a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
+++ b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
@@ -27,7 +27,8 @@ function main()
   η_coeff = 2
   α_coeff = 4
   vf = 0.4
-  path = dirname(dirname(@__DIR__))*"/results/nonlinear_thermal_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/nonlinear_thermal_compliance_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -88,11 +89,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it, uh, φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/Serial/thermal_compliance_ALM.jl b/scripts/Serial/thermal_compliance_ALM.jl
index 171d8d49..425a63d5 100644
--- a/scripts/Serial/thermal_compliance_ALM.jl
+++ b/scripts/Serial/thermal_compliance_ALM.jl
@@ -26,7 +26,8 @@ function main()
   vf = 0.4
   η_coeff = 2
   α_coeff = 4
-  path = dirname(dirname(@__DIR__))*"/results/thermal_compliance_ALM"
+  path = dirname(dirname(@__DIR__))*"/results/thermal_compliance_ALM/"
+  iter_mod = 10
   mkdir(path)
 
   ## FE Setup
@@ -86,11 +87,12 @@ function main()
   optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
+    iszero(it % iter_mod) && writevtk(Ω,path*"out$it",cellfields=data)
     write_history(path*"/history.txt",optimiser.history)
   end
   it = get_history(optimiser).niter; uh = get_state(pcfs)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 main()
\ No newline at end of file
diff --git a/scripts/_dev/3d_hyperelastic_compliance_ALM.jl b/scripts/_dev/3d_hyperelastic_compliance_ALM.jl
index 8066e3e6..2565536b 100644
--- a/scripts/_dev/3d_hyperelastic_compliance_ALM.jl
+++ b/scripts/_dev/3d_hyperelastic_compliance_ALM.jl
@@ -73,7 +73,7 @@ function main(mesh_partition,distribute,el_size,δₓ)
   ## Optimiser
   u = LevelSetTopOpt.forward_solve!(state_map,φh)
   uh = FEFunction(U,u)
-  write_vtk(Ω,path*"/struc_$δₓ",0,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"/struc_$δₓ",cellfields=["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 function main_alt(mesh_partition,distribute,el_size,gz)
@@ -173,7 +173,7 @@ function main_alt(mesh_partition,distribute,el_size,gz)
   ## Optimiser
   u = LevelSetTopOpt.forward_solve!(state_map,φh)
   uh = FEFunction(U,u)
-  write_vtk(Ω,path*"/struc_neohook_$gz",0,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
+  writevtk(Ω,path*"/struc_neohook_$gz",cellfields=["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
 end
 
 with_mpi() do distribute
diff --git a/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl b/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl
index aa2d019e..4ab969a6 100644
--- a/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl
+++ b/scripts/_dev/inv_hom_block_assem_testing_3D_MPI.jl
@@ -136,13 +136,14 @@ function main(mesh_partition,distribute,el_size,diag_assem::Bool)
     λi = optimiser.λ; Λi = optimiser.Λ
     print_history(it,["J"=>Ji,"C"=>Ci,"L"=>Li,"λ"=>λi,"Λ"=>Λi];ranks=ranks)
     write_history(history,path*"/history.csv";ranks=ranks)
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]
+    iszero(it % 10) && writevtk(Ω,path*"out$it",cellfields=data)
   end
   it,Ji,Ci,Li = last(optimiser.history)
   λi = optimiser.λ; Λi = optimiser.Λ
   print_history(it,["J"=>Ji,"C"=>Ci,"L"=>Li,"λ"=>λi,"Λ"=>Λi];ranks=ranks)
   write_history(optimiser.history,path*"/history.csv";ranks=ranks)
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))];iter_mod=1)
+  writevtk(Ω,path*"/out$it",cellfields=["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
 end
 
 # RUN: mpiexecjl --project=. -n 8 julia ./scripts/inv_hom_block_assem_testing_3D_MPI.jl
diff --git a/scripts/_dev/inv_hom_block_assem_testing_full_script.jl b/scripts/_dev/inv_hom_block_assem_testing_full_script.jl
index 75018ee6..f3c31fc2 100644
--- a/scripts/_dev/inv_hom_block_assem_testing_full_script.jl
+++ b/scripts/_dev/inv_hom_block_assem_testing_full_script.jl
@@ -115,13 +115,14 @@ function main(diag_block)
     λi = optimiser.λ; Λi = optimiser.Λ
     print_history(it,["J"=>Ji,"C"=>Ci,"L"=>Li,"λ"=>λi,"Λ"=>Λi])
     write_history(history,path*"/history.csv")
-    write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
+    data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|∇(φ)|"=>(norm ∘ ∇(φh))]
+    iszero(it % 10) && writevtk(Ω,path*"out$it",cellfields=data)
   end
   it,Ji,Ci,Li = last(optimiser.history)
   λi = optimiser.λ; Λi = optimiser.Λ
   print_history(it,["J"=>Ji,"C"=>Ci,"L"=>Li,"λ"=>λi,"Λ"=>Λi])
   write_history(optimiser.history,path*"/history.csv")
-  write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))];iter_mod=1)
+  writevtk(Ω,path*"/out$it",cellfields=["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh))])
 end
 
 using LinearAlgebra
diff --git a/src/LevelSetTopOpt.jl b/src/LevelSetTopOpt.jl
index 39917e95..7b6ebc4d 100644
--- a/src/LevelSetTopOpt.jl
+++ b/src/LevelSetTopOpt.jl
@@ -51,7 +51,6 @@ export update_labels!
 export initial_lsf
 export get_el_Δ
 export isotropic_elast_tensor
-export write_vtk
 
 include("Advection.jl")
 export AdvectionStencil
diff --git a/src/Optimisers/AugmentedLagrangian.jl b/src/Optimisers/AugmentedLagrangian.jl
index fb9d74fe..b2380686 100644
--- a/src/Optimisers/AugmentedLagrangian.jl
+++ b/src/Optimisers/AugmentedLagrangian.jl
@@ -178,6 +178,10 @@ function Base.iterate(m::AugmentedLagrangian,state)
   update_mod, reinit_mod, ζ, Λ_max, γ_reinit, os_γ_mult = params.update_mod, 
     params.reinit_mod, params.ζ, params.Λ_max, params.γ_reinit, params.os_γ_mult
 
+  ## Periodicially call GC
+  iszero(it % 50) && GC.gc();
+
+  ## Check stopping criteria
   if finished(m)
     return nothing
   end
diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl
index c64a9d7f..9cb75952 100644
--- a/src/Optimisers/HilbertianProjection.jl
+++ b/src/Optimisers/HilbertianProjection.jl
@@ -321,6 +321,10 @@ function Base.iterate(m::HilbertianProjection,state)
   it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it = state
   history, params = m.history, m.params
 
+  ## Periodicially call GC
+  iszero(it % 50) && GC.gc();
+
+  ## Check stopping criteria
   if finished(m)
     return nothing
   end
diff --git a/src/Utilities.jl b/src/Utilities.jl
index 82ba9f0f..aede0bd2 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -212,29 +212,4 @@ function isotropic_elast_tensor(D::Int,E::Number,v::Number)
   else
     @notimplemented
   end
-end
-
-"""
-    write_vtk(Ω,path,it,entries::Vector{<:Pair};iter_mod=10)
-
-Write a VTK file to `path`. This functions in a similar way to
-Gridap's `writevtk` function except we 
-
-!!! note
-    This may be removed in future and replaced by
-
-    ```(isone(it) || iszero(it % iter_mod)) && writevtk(Ω,path,cellfields=entries)```
-
-    or 
-    
-    ```iszero(it-1 % iter_mod) && writevtk(Ω,path,cellfields=entries)```
-"""
-function write_vtk(Ω,path,it,entries::Vector{<:Pair};iter_mod=10) # TODO: Rename to writevtk?
-  if isone(it) || iszero(it % iter_mod) 
-    writevtk(Ω,path,cellfields=entries)
-  end
-
-  if iszero(it % 10*iter_mod)
-    GC.gc() # TODO: Test in Julia 1.10.0 and check if fixed.
-  end 
-end
+end
\ No newline at end of file