diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 561fd70c74bedf4b3c1c26c2e9115a7d3044706c..5baeae1f53827be4226455eda278c6d5580f39d1 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -474,9 +474,8 @@ void TimeLoop::initialize()
 
     // Output initial conditions
     {
-        const bool output_initial_condition = true;
-        outputSolutions(output_initial_condition, 0, _start_time,
-                        &Output::doOutput);
+        preOutputInitialConditions(_start_time);
+        outputSolutions(0, _start_time, &Output::doOutput);
     }
 
     auto const time_step_constraints = generateOutputTimeStepConstraints(
@@ -530,9 +529,7 @@ bool TimeLoop::calculateNextTimeStep()
 
     if (!_last_step_rejected)
     {
-        const bool output_initial_condition = false;
-        outputSolutions(output_initial_condition, timesteps, current_time,
-                        &Output::doOutput);
+        outputSolutions(timesteps, current_time, &Output::doOutput);
     }
 
     if (std::abs(_current_time - _end_time) <
@@ -564,9 +561,7 @@ void TimeLoop::outputLastTimeStep() const
     // output last time step
     if (successful_time_step)
     {
-        const bool output_initial_condition = false;
-        outputSolutions(output_initial_condition,
-                        _accepted_steps + _rejected_steps, _current_time,
+        outputSolutions(_accepted_steps + _rejected_steps, _current_time,
                         &Output::doOutputLastTimestep);
     }
 }
@@ -784,22 +779,16 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
     return nonlinear_solver_status;
 }
 
-void TimeLoop::outputSolutions(bool const output_initial_condition) const
+void TimeLoop::outputSolutions() const
 {
     const std::size_t timesteps = _accepted_steps + 1;
-    outputSolutions(output_initial_condition, timesteps, _current_time,
-                    &Output::doOutput);
+    outputSolutions(timesteps, _current_time, &Output::doOutput);
 }
 
 template <typename OutputClassMember>
-void TimeLoop::outputSolutions(bool const output_initial_condition,
-                               unsigned timestep, const double t,
+void TimeLoop::outputSolutions(unsigned timestep, const double t,
                                OutputClassMember output_class_member) const
 {
-    // All _per_process_data share the first process.
-    bool const is_staggered_coupling =
-        !isMonolithicProcess(*_per_process_data[0]);
-
     for (auto const& process_data : _per_process_data)
     {
         // If nonlinear solver diverged, the solution has already been
@@ -810,45 +799,8 @@ void TimeLoop::outputSolutions(bool const output_initial_condition,
         }
 
         auto const process_id = process_data->process_id;
-        auto& pcs = process_data->process;
+        auto const& pcs = process_data->process;
 
-        if (!is_staggered_coupling && output_initial_condition)
-        {
-            // dummy value to handle the time derivative terms more or less
-            // correctly, i.e. to ignore them.
-            double const dt = 1;
-            process_data->time_disc->nextTimestep(t, dt);
-
-            pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
-            // Update secondary variables, which might be uninitialized, before
-            // output.
-            pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
-                                         *_process_solutions_prev[process_id],
-                                         process_id);
-        }
-        if (is_staggered_coupling && output_initial_condition)
-        {
-            CoupledSolutionsForStaggeredScheme coupled_solutions(
-                _process_solutions);
-
-            process_data->process.setCoupledSolutionsForStaggeredScheme(
-                &coupled_solutions);
-            process_data->process
-                .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
-                    process_id);
-
-            // dummy value to handle the time derivative terms more or less
-            // correctly, i.e. to ignore them.
-            double const dt = 1;
-            process_data->time_disc->nextTimestep(t, dt);
-
-            pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
-            // Update secondary variables, which might be uninitialized, before
-            // output.
-            pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
-                                         *_process_solutions_prev[process_id],
-                                         process_id);
-        }
         for (auto const& output_object : _outputs)
         {
             (output_object.*output_class_member)(
@@ -902,5 +854,63 @@ double TimeLoop::computeRelativeSolutionChangeFromPreviousTimestep(
     const double solution_error =
         NumLib::computeRelativeChangeFromPreviousTimestep(x, x_prev, norm_type);
     return solution_error;
-};
+}
+
+void TimeLoop::preOutputInitialConditions(const double t) const
+{
+    // All _per_process_data share the first process.
+    bool const is_staggered_coupling =
+        !isMonolithicProcess(*_per_process_data[0]);
+
+    for (auto const& process_data : _per_process_data)
+    {
+        // If nonlinear solver diverged, the solution has already been
+        // saved.
+        if (!process_data->nonlinear_solver_status.error_norms_met)
+        {
+            continue;
+        }
+
+        auto const process_id = process_data->process_id;
+        auto& pcs = process_data->process;
+
+        if (!is_staggered_coupling)
+        {
+            // dummy value to handle the time derivative terms more or less
+            // correctly, i.e. to ignore them.
+            double const dt = 1;
+            process_data->time_disc->nextTimestep(t, dt);
+
+            pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
+            // Update secondary variables, which might be uninitialized, before
+            // output.
+            pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
+                                         *_process_solutions_prev[process_id],
+                                         process_id);
+        }
+        else
+        {
+            CoupledSolutionsForStaggeredScheme coupled_solutions(
+                _process_solutions);
+
+            process_data->process.setCoupledSolutionsForStaggeredScheme(
+                &coupled_solutions);
+            process_data->process
+                .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
+                    process_id);
+
+            // dummy value to handle the time derivative terms more or less
+            // correctly, i.e. to ignore them.
+            double const dt = 1;
+            process_data->time_disc->nextTimestep(t, dt);
+
+            pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
+            // Update secondary variables, which might be uninitialized, before
+            // output.
+            pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
+                                         *_process_solutions_prev[process_id],
+                                         process_id);
+        }
+    }
+}
 }  // namespace ProcessLib
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index 62d09f3b6a3f682e45838e812508f0168ee1befa..1e83f163116b8dbd5b02f101fb3f66b7c2f849c2 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -62,7 +62,7 @@ public:
     double endTime() const { return _end_time; }
     double currentTime() const { return _current_time; }
     bool successful_time_step = false;
-    void outputSolutions(bool const output_initial_condition) const;
+    void outputSolutions() const;
 
 private:
     bool preTsNonlinearSolvePostTs(double const t, double const dt,
@@ -124,7 +124,7 @@ private:
             time_step_constraints);
 
     template <typename OutputClassMember>
-    void outputSolutions(bool const output_initial_condition, unsigned timestep,
+    void outputSolutions(unsigned timestep,
                          const double t,
                          OutputClassMember output_class_member) const;
 
@@ -133,6 +133,7 @@ private:
     generateOutputTimeStepConstraints(std::vector<double>&& fixed_times) const;
     double computeRelativeSolutionChangeFromPreviousTimestep(
         double const t, std::size_t process_index) const;
+    void preOutputInitialConditions(const double t) const;
     std::vector<GlobalVector*> _process_solutions;
     std::vector<GlobalVector*> _process_solutions_prev;
     std::vector<Output> _outputs;