diff --git a/src/program_state.hh b/src/program_state.hh
index f05e193cd738ce8f772995c56f9d2a02c793252b..92aeb766048326a2e5f0632b25207cf4d650dd1d 100644
--- a/src/program_state.hh
+++ b/src/program_state.hh
@@ -68,7 +68,7 @@ template <class Vector, class ScalarVector> class ProgramState {
       solver.solve();
     };
 
-    timeStep = 1;
+    timeStep = 0;
     relativeTime = 0.0;
     relativeTau = 1e-6;
 
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 3942ef1fe73e575fe310e682f5151a4584a3d723..bd7b8182f9493797e0ed3af45da1e33a9b29e847 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -178,12 +178,6 @@ int main(int argc, char *argv[]) {
                    functions.get("velocityDirichletCondition"),
                &neumannFunction = functions.get("neumannCondition");
 
-    std::fstream timeStepWriter("timeSteps", std::fstream::out);
-    timeStepWriter << std::fixed << std::setprecision(10);
-    auto const reportTimeStep = [&](double _relativeTime, double _relativeTau) {
-      timeStepWriter << _relativeTime << " " << _relativeTau << std::endl;
-    };
-
     MyAssembler const myAssembler(leafView);
 
     MyBody<dims> const body(parset);
@@ -237,6 +231,7 @@ int main(int argc, char *argv[]) {
       }
     }
 
+    // Set up writers
     FrictionWriter<ScalarVector, Vector> frictionWriter(
         vertexCoordinates, frictionalNodes, "friction",
         MyGeometry::horizontalProjection);
@@ -251,37 +246,43 @@ int main(int argc, char *argv[]) {
     SpecialWriter<GridView, dims> specialDisplacementWriter(
         "specialDisplacements", leafView);
 
-    auto const report = [&](Vector const &_u, Vector const &_v,
-                            ScalarVector const &_alpha) {
-      horizontalSurfaceWriter.writeKinetics(_u, _v);
-      verticalSurfaceWriter.writeKinetics(_u, _v);
+    std::fstream timeStepWriter("timeSteps", std::fstream::out);
+    timeStepWriter << std::fixed << std::setprecision(10);
+
+    MyVTKWriter<typename MyAssembler::VertexBasis,
+                typename MyAssembler::CellBasis> const
+        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
+    auto const report = [&]() {
+      if (programState.timeStep != 0)
+        timeStepWriter << programState.relativeTime << " "
+                       << programState.relativeTau << std::endl;
+
+      horizontalSurfaceWriter.writeKinetics(programState.u, programState.v);
+      verticalSurfaceWriter.writeKinetics(programState.u, programState.v);
 
       ScalarVector c;
-      myGlobalFriction->coefficientOfFriction(_v, c);
-      frictionWriter.writeKinetics(_u, _v);
-      frictionWriter.writeOther(c, _alpha);
+      myGlobalFriction->coefficientOfFriction(programState.v, c);
+      frictionWriter.writeKinetics(programState.u, programState.v);
+      frictionWriter.writeOther(c, programState.alpha);
 
       BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
-          myAssembler.vertexBasis, _v);
+          myAssembler.vertexBasis, programState.v);
       BasisGridFunction<typename MyAssembler::VertexBasis, Vector> displacement(
-          myAssembler.vertexBasis, _u);
+          myAssembler.vertexBasis, programState.u);
       specialVelocityWriter.write(velocity);
       specialDisplacementWriter.write(displacement);
-    };
-    report(programState.u, programState.v, programState.alpha);
-
-    MyVTKWriter<typename MyAssembler::VertexBasis,
-                typename MyAssembler::CellBasis> const
-        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
-    if (parset.get<bool>("io.writeVTK")) {
-      ScalarVector stress;
-      myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                         body.getPoissonRatio(), programState.u,
-                                         stress);
-      vtkWriter.write(0, programState.u, programState.v, programState.alpha,
-                      stress);
-    }
+      if (parset.get<bool>("io.writeVTK")) {
+        ScalarVector stress;
+        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+                                           body.getPoissonRatio(),
+                                           programState.u, stress);
+        vtkWriter.write(programState.timeStep, programState.u, programState.v,
+                        programState.alpha, stress);
+      }
+    };
+    report();
 
     // Set up TNNMG solver
     using NonlinearFactory = SolverFactory<
@@ -322,25 +323,17 @@ int main(int argc, char *argv[]) {
                             programState.relativeTime, programState.relativeTau,
                             computeExternalForces, stateEnergyNorm, mustRefine);
     while (!adaptiveTimeStepper.reachedEnd()) {
+      programState.timeStep++;
+
       adaptiveTimeStepper.advance();
+
       programState.relativeTime = adaptiveTimeStepper.getRelativeTime();
       programState.relativeTau = adaptiveTimeStepper.getRelativeTau();
-      reportTimeStep(programState.relativeTime, programState.relativeTau);
-
       current.second->extractDisplacement(programState.u);
       current.second->extractVelocity(programState.v);
       current.first->extractAlpha(programState.alpha);
 
-      report(programState.u, programState.v, programState.alpha);
-      if (parset.get<bool>("io.writeVTK")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(),
-                                           programState.u, stress);
-        vtkWriter.write(programState.timeStep, programState.u, programState.v,
-                        programState.alpha, stress);
-      }
-      programState.timeStep++;
+      report();
     }
     timeStepWriter.close();
     Python::stop();