diff --git a/examples/dem.py b/examples/dem.py
index 4c282f4ad9bc78263eb8fdc876348fff6056ff79..33a4e366756b7419884a95cd834b856b8c9b3900 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -162,7 +162,7 @@ psim.read_particle_data(
     pairs.halfspace())
 
 psim.build_neighbor_lists(linkedCellWidth + skin)
-psim.vtk_output(f"output/dem_{target}")
+psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
 psim.compute(linear_spring_dashpot, linkedCellWidth + skin, symbols={'dt': dt_SI})
 psim.compute(euler, symbols={'dt': dt_SI})
 psim.compute(gravity, symbols={'densityParticle_SI': densityParticle_SI,
diff --git a/runtime/vtk.hpp b/runtime/vtk.hpp
index e2041a787dd0468e4a37995721f1e731ebc52c05..ddfb659d915cdb299d4fbffa98567e13ae40eba2 100644
--- a/runtime/vtk.hpp
+++ b/runtime/vtk.hpp
@@ -8,13 +8,17 @@
 
 namespace pairs {
 
-void vtk_write_data(PairsSimulation *ps, const char *filename, int start, int end, int timestep) {
+void vtk_write_data(PairsSimulation *ps, const char *filename, int start, int end, int timestep, int frequency) {
     std::string output_filename(filename);
     auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
     auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
     const int n = end - start;
     std::ostringstream filename_oss;
 
+    if(frequency != 0 && timestep % frequency != 0) {
+        return;
+    }
+
     filename_oss << filename << "_";
     if(ps->getDomainPartitioner()->getWorldSize() > 1) {
         filename_oss << "r" << ps->getDomainPartitioner()->getRank() << "_";
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index edc321f51188be13435949fb754de0734076c9dd..6315036dd311555cbb5ddcc00f44ebc3059c7fa5 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -69,6 +69,7 @@ class Simulation:
         self.expr_id = 0
         self.iter_id = 0
         self.vtk_file = None
+        self.vtk_frequency = 0
         self._target = None
         self._dom_part = DimensionRanges(self)
 
@@ -247,8 +248,9 @@ class Simulation:
         else:
             self.nested_count += 1
 
-    def vtk_output(self, filename):
+    def vtk_output(self, filename, frequency=0):
         self.vtk_file = filename
+        self.vtk_frequency = frequency
 
     def target(self, target):
         self._target = target
@@ -283,13 +285,13 @@ class Simulation:
 
         timestep = Timestep(self, self.ntimesteps, timestep_procedures)
         self.enter(timestep.block)
-        timestep.add(VTKWrite(self, self.vtk_file, timestep.timestep() + 1))
+        timestep.add(VTKWrite(self, self.vtk_file, timestep.timestep() + 1, self.vtk_frequency))
         self.leave()
 
         body = Block.from_list(self, [
             self.setups,
             BuildCellListsStencil(self, self.cell_lists),
-            VTKWrite(self, self.vtk_file, 0),
+            VTKWrite(self, self.vtk_file, 0, self.vtk_frequency),
             timestep.as_block()
         ])
 
diff --git a/src/pairs/sim/vtk.py b/src/pairs/sim/vtk.py
index 806afc82d0705a7e813407139c677acae24e35e2..68f84a1fd0257bf769929449ee8219a9b5930137 100644
--- a/src/pairs/sim/vtk.py
+++ b/src/pairs/sim/vtk.py
@@ -6,15 +6,16 @@ from pairs.sim.lowerable import Lowerable
 
 
 class VTKWrite(Lowerable):
-    def __init__(self, sim, filename, timestep):
+    def __init__(self, sim, filename, timestep, frequency):
         super().__init__(sim)
         self.filename = filename
         self.timestep = Lit.cvt(sim, timestep)
+        self.frequency = frequency
 
     @pairs_inline
     def lower(self):
         nlocal = self.sim.nlocal
         nghost = self.sim.nghost
         nall = nlocal + nghost
-        Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_local", 0, nlocal, self.timestep])
-        Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_ghost", nlocal, nall, self.timestep])
+        Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_local", 0, nlocal, self.timestep, self.frequency])
+        Call_Void(self.sim, "pairs::vtk_write_data", [self.filename + "_ghost", nlocal, nall, self.timestep, self.frequency])