diff --git a/examples/adjoints/adjoints_ex5/adjoints_ex5.C b/examples/adjoints/adjoints_ex5/adjoints_ex5.C
index f50d3de3787..4831307e4c4 100644
--- a/examples/adjoints/adjoints_ex5/adjoints_ex5.C
+++ b/examples/adjoints/adjoints_ex5/adjoints_ex5.C
@@ -642,9 +642,6 @@ int main (int argc, char ** argv)
       // A SensitivityData object
       SensitivityData sensitivities(qois, system, system.get_parameter_vector());
 
-      // Accumulator for the time integrated total sensitivity
-      Number total_sensitivity = 0.0;
-
       // Retrieve the primal and adjoint solutions at the current timestep
       system.time_solver->retrieve_timestep();
 
@@ -673,6 +670,26 @@ int main (int argc, char ** argv)
                    << std::endl
                    << std::endl;
 
+      ParameterVector & parameter_vector = system.get_parameter_vector();
+
+      auto integrate_adjoint_sensitivity = [&qois, &parameter_vector, &sensitivities](Real time_quadrature_weight, System & system)
+      {
+        SensitivityData sensitivity_contributions(qois, system, parameter_vector);
+
+        system.adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivity_contributions);
+
+        // Remove the sensitivity rhs vector from system since we did not write it to file and it cannot be retrieved
+        system.remove_vector("sensitivity_rhs0");
+
+        // Get the contributions for each sensitivity from this timestep
+        for(unsigned int i = 0; i != qois.size(system); i++)
+          for(unsigned int j = 0; j != parameter_vector.size(); j++)
+            sensitivities[i][j] += ( sensitivity_contributions[i][j] )*(time_quadrature_weight);
+      };
+
+      std::vector<std::function<void(Real, System &)>> sensitivity_operations;
+      sensitivity_operations.push_back(integrate_adjoint_sensitivity);
+
       // Now we begin the timestep loop to compute the time-accurate
       // adjoint sensitivities
       for (unsigned int t_step=param.initial_timestep;
@@ -680,7 +697,7 @@ int main (int argc, char ** argv)
         {
           // Call the postprocess function which we have overloaded to compute
           // accumulate the perturbed residuals
-          system.time_solver->integrate_adjoint_sensitivity(qois, system.get_parameter_vector(), sensitivities);
+          system.time_solver->advance_postprocessing_timestep(sensitivity_operations);
 
           // A pretty update message
           libMesh::out << "Retrieved, "
@@ -706,14 +723,11 @@ int main (int argc, char ** argv)
                    << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
                    << std::endl
                    << std::endl;
-
-          // Add the contribution of the current timestep to the total sensitivity
-          total_sensitivity += sensitivities[0][0];
         }
 
       // Print it out
       libMesh::out << "Sensitivity of QoI 0 w.r.t parameter 0 is: "
-                   << total_sensitivity
+                   << sensitivities[0][0]
                    << std::endl;
 
       // Hard coded test to ensure that the actual numbers we are
@@ -725,12 +739,12 @@ int main (int argc, char ** argv)
         libmesh_error_msg_if(std::abs(system.time - (1.0089)) >= 2.e-4,
                              "Mismatch in end time reached by adaptive timestepper!");
 
-        libmesh_error_msg_if(std::abs(total_sensitivity - 4.87767) >= 3.e-3,
+        libmesh_error_msg_if(std::abs(sensitivities[0][0] - 4.87767) >= 3.e-3,
                              "Mismatch in sensitivity gold value!");
       }
       else
       {
-        libmesh_error_msg_if(std::abs(total_sensitivity - 4.83551) >= 2.e-4,
+        libmesh_error_msg_if(std::abs(sensitivities[0][0] - 4.83551) >= 2.e-4,
                              "Mismatch in sensitivity gold value!");
       }
 #ifdef NDEBUG
diff --git a/examples/adjoints/adjoints_ex7/adjoints_ex7.C b/examples/adjoints/adjoints_ex7/adjoints_ex7.C
index 47cc6fa2c8c..a987cdd45a9 100644
--- a/examples/adjoints/adjoints_ex7/adjoints_ex7.C
+++ b/examples/adjoints/adjoints_ex7/adjoints_ex7.C
@@ -460,9 +460,6 @@ int main (int argc, char ** argv)
   mesh.print_info();
   equation_systems.print_info();
 
-  // Accumulated integrated QoIs
-  Number QoI_1_accumulated = 0.0;
-
   // In optimized mode we catch any solver errors, so that we can
   // write the proper footers before closing.  In debug mode we just
   // let the exception throw so that gdb can grab it.
@@ -511,18 +508,39 @@ int main (int argc, char ** argv)
 
       system.time_solver->retrieve_timestep();
 
-      QoI_1_accumulated = 0.0;
+      std::vector<Number> QoI_accumulated(system.n_qois(), 0.0);
+
+      // Create a lambda to integrate a qoi across a timestep
+      auto integrate_qoi_timestep = [&QoI_accumulated] (Real time_quadrature_weight, System & system)
+      {
+        // Zero out the system.qoi vector
+        for (auto j : make_range(system.n_qois()))
+        {
+         (system.qoi)[j] = 0.0;
+        }
+
+        system.assemble_qoi();
+
+        for (auto j : make_range(system.n_qois()))
+        {
+          QoI_accumulated[j] += time_quadrature_weight*((system.qoi)[j]);
+        }
+
+        return;
+      };
+
+      std::vector<std::function<void(Real, System &)>> integration_operations;
+
+      integration_operations.push_back(integrate_qoi_timestep);
 
       for (unsigned int t_step=param.initial_timestep;
            t_step != param.initial_timestep + param.n_timesteps; ++t_step)
         {
-          system.time_solver->integrate_qoi_timestep();
-
-          QoI_1_accumulated += (system.qoi)[1];
+          system.time_solver->advance_postprocessing_timestep(integration_operations);
         }
 
-      std::cout<< "The computed QoI 0 is " << std::setprecision(17) << (system.qoi)[0] << std::endl;
-      std::cout<< "The computed QoI 1 is " << std::setprecision(17) << QoI_1_accumulated << std::endl;
+      std::cout<< "The computed QoI 0 is " << std::setprecision(17) << QoI_accumulated[0] << std::endl;
+      std::cout<< "The computed QoI 1 is " << std::setprecision(17) << QoI_accumulated[1] << std::endl;
 
       ///////////////// Now for the Adjoint Solution //////////////////////////////////////
 
@@ -646,15 +664,6 @@ int main (int argc, char ** argv)
           // unnecessarily in the error estimator
           system.set_adjoint_already_solved(true);
 
-          libMesh::out << "Saving adjoint and retrieving primal solutions at time t=" << system.time - system.deltat << std::endl;
-
-          // The adjoint_advance_timestep function calls the retrieve and store
-          // function of the memory_solution_history class via the
-          // memory_solution_history object we declared earlier.  The
-          // retrieve function sets the system primal vectors to their
-          // values at the current time instance.
-          system.time_solver->adjoint_advance_timestep();
-
           libMesh::out << "|Z0("
                        << system.time
                        << ")|= "
@@ -669,6 +678,15 @@ int main (int argc, char ** argv)
                        << std::endl
                        << std::endl;
 
+          libMesh::out << "Saving adjoint and retrieving primal solutions at time t=" << system.time << std::endl;
+
+          // The adjoint_advance_timestep function calls the retrieve and store
+          // function of the memory_solution_history class via the
+          // memory_solution_history object we declared earlier.  The
+          // retrieve function sets the system primal vectors to their
+          // values at the current time instance.
+          system.time_solver->adjoint_advance_timestep();
+
           // Get a pointer to the primal solution vector
           primal_solution = *system.solution;
 
@@ -704,8 +722,8 @@ int main (int argc, char ** argv)
       // as we did for the adjoint solve. This object will now define the residual for the dual weighted error evaluation at
       // each timestep.
       // For adjoint consistency, it is important that we maintain the same integration method for the error estimation integral
-      // as we did for the primal and QoI integration. This is ensured within the library, but is currently ONLY supported for
-      // the Backward-Euler (theta = 0.5) time integration method.
+      // as we did for the primal and QoI integration. This is ensured within the library, but is currently only supported for
+      // the Backward-Euler (theta = 1.0) time integration method.
       QoISet qois;
 
       std::vector<unsigned int> qoi_indices;
@@ -719,55 +737,94 @@ int main (int argc, char ** argv)
       auto adjoint_refinement_error_estimator =
         build_adjoint_refinement_error_estimator(qois, dynamic_cast<FEMPhysics *>(sigma_physics), param);
 
-      // Error Vector to be filled by the error estimator at each timestep
+      std::vector<Number> time_integrated_QoI_error(system.n_qois(), 0.0);
+      system.qoi_error_estimates.resize(system.n_qois());
+      std::vector< std::unique_ptr<NumericVector<Number>> > old_adjoints;
       ErrorVector QoI_elementwise_error;
-      std::vector<Number> QoI_spatially_integrated_error(system.n_qois());
 
-      // Error Vector to hold the total accumulated error
-      ErrorVector accumulated_QoI_elementwise_error;
-      std::vector<Number> accumulated_QoI_spatially_integrated_error(system.n_qois());
+      old_adjoints.resize(system.n_qois());
 
-      // Retrieve the primal and adjoint solutions at the current timestep
-      system.time_solver->retrieve_timestep();
+      // Initialize old adjoints
+      for(auto j : make_range(system.n_qois()))
+      {
+       old_adjoints[j] = system.get_adjoint_solution(j).clone();
+      }
 
-      // A pretty update message
-      libMesh::out << "Retrieved, "
-               << "time = "
-               << system.time
-               << std::endl;
+      auto integrate_adjoint_refinement_error_estimate = [&QoI_elementwise_error, &time_integrated_QoI_error, &old_adjoints, &adjoint_refinement_error_estimator] (Real time_quadrature_weight, System & system)
+      {
+        // We will need to move the system.time around to ensure that residuals
+        // are built with the right deltat and the right time.
+        Real next_step_deltat;
 
-      libMesh::out << "|U("
-                   << system.time
-                   << ")|= "
-                   << system.calculate_norm(*system.solution, 0, H1)
-                   << std::endl;
+        // Reset the elementwise QoI error vector
+        QoI_elementwise_error.clear();
 
-      libMesh::out << "|U_old("
-                   << system.time
-                   << ")|= "
-                   << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
-                   << std::endl;
+        // Zero out the system.qoi vector
+        for (auto j : make_range(system.n_qois()))
+        {
+         (system.qoi_error_estimates)[j] = 0.0;
+        }
 
-      libMesh::out << "|Z0("
-                   << system.time
-                   << ")|= "
-                   << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
-                   << std::endl
-                   << std::endl;
+        // For evaluating the residual, we need to use the deltat that was used
+        // to get us to this solution, so we save the current deltat as next_step_deltat
+        // and set _system.deltat to the last completed deltat.
+        next_step_deltat = dynamic_cast<HeatSystem &>(system).deltat;
+        dynamic_cast<HeatSystem &>(system).deltat = dynamic_cast<HeatSystem &>(system).get_time_solver().get_last_deltat();
 
-      libMesh::out << "|Z1("
-                   << system.time
-                   << ")|= "
-                   << system.calculate_norm(system.get_adjoint_solution(1), 0, H1)
-                   << std::endl
-                   << std::endl;
+        // The adjoint error estimate expression for a backward Euler
+        // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
+        for (auto j : make_range(system.n_qois()))
+        {
+          // Swap for residual weighting
+          system.get_adjoint_solution(j).swap(*old_adjoints[j]);
+        }
+
+        system.update();
+
+        // The residual has to be evaluated at the last time
+        system.time = system.time - dynamic_cast<HeatSystem &>(system).get_time_solver().get_last_deltat();
+
+        adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
+
+        for (auto j : make_range(system.n_qois()))
+        {
+          (system.qoi_error_estimates)[j] += adjoint_refinement_error_estimator->get_global_QoI_error_estimate(j);
+          //std::cout<<"Error est: "<<(system.qoi_error_estimates)[j]<<std::endl;
+          time_integrated_QoI_error[j] += time_quadrature_weight*((system.qoi_error_estimates)[j]);
+        }
+
+        // Shift the time back
+        system.time = system.time + dynamic_cast<HeatSystem &>(system).get_time_solver().get_last_deltat();
+
+        // Swap back the current and old adjoints
+        for (auto j : make_range(system.n_qois()))
+        {
+          system.get_adjoint_solution(j).swap(*old_adjoints[j]);
+        }
+
+        system.update();
+
+        // Set the system deltat back to what it should be to march to the next time
+        dynamic_cast<HeatSystem &>(system).deltat = next_step_deltat;
+
+        // The adjoint error estimate expression for a backwards Euler
+        // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
+        for (auto j : make_range(system.n_qois()))
+        {
+          old_adjoints[j] = system.get_adjoint_solution(j).clone();
+        }
+      };
+
+      std::vector<std::function<void(Real, System &)>> error_estimation_operations;
+
+      error_estimation_operations.push_back(integrate_adjoint_refinement_error_estimate);
 
       // Now we begin the timestep loop to compute the time-accurate
       // adjoint sensitivities
       for (unsigned int t_step=param.initial_timestep;
            t_step != param.initial_timestep + param.n_timesteps; ++t_step)
         {
-          system.time_solver->integrate_adjoint_refinement_error_estimate(*adjoint_refinement_error_estimator, QoI_elementwise_error);
+          system.time_solver->advance_postprocessing_timestep(error_estimation_operations);
 
           // A pretty update message
           libMesh::out << "Retrieved, "
@@ -801,23 +858,10 @@ int main (int argc, char ** argv)
                    << std::endl
                    << std::endl;
 
-          // Error contribution from this timestep
-          for(unsigned int i = 0; i < QoI_elementwise_error.size(); i++)
-            accumulated_QoI_elementwise_error[i] += QoI_elementwise_error[i];
-
-          // QoI wise error contribution from this timestep
-          for (auto j : make_range(system.n_qois()))
-          {
-            // Skip this QoI if not in the QoI Set
-            if ((adjoint_refinement_error_estimator->qoi_set()).has_index(j))
-            {
-              accumulated_QoI_spatially_integrated_error[j] += (system.qoi_error_estimates)[j];
-            }
-          }
         }
 
-      std::cout<<"Time integrated error estimate for QoI 0: "<<std::setprecision(17)<<accumulated_QoI_spatially_integrated_error[0]<<std::endl;
-      std::cout<<"Time integrated error estimate for QoI 1: "<<std::setprecision(17)<<accumulated_QoI_spatially_integrated_error[1]<<std::endl;
+      std::cout<<"Time integrated error estimate for QoI 0: "<<std::setprecision(17)<<time_integrated_QoI_error[0]<<std::endl;
+      std::cout<<"Time integrated error estimate for QoI 1: "<<std::setprecision(17)<<time_integrated_QoI_error[1]<<std::endl;
 
       // Hard coded test to ensure that the actual numbers we are
       // getting are what they should be
@@ -828,18 +872,18 @@ int main (int argc, char ** argv)
         libmesh_error_msg_if(std::abs(system.time - (1.7548735069535084)) >= 2.e-4,
                              "Mismatch in end time reached by adaptive timestepper!");
 
-        libmesh_error_msg_if(std::abs(accumulated_QoI_spatially_integrated_error[0] - (-0.75139371165754232)) >= 2.e-4,
+        libmesh_error_msg_if(std::abs(time_integrated_QoI_error[0] - (-0.75139371165754232)) >= 2.e-4,
                              "Error Estimator identity not satisfied!");
 
-        libmesh_error_msg_if(std::abs(accumulated_QoI_spatially_integrated_error[1] - (-0.70905030091469889)) >= 2.e-4,
+        libmesh_error_msg_if(std::abs(time_integrated_QoI_error[1] - (-0.70905030091469889)) >= 2.e-4,
                              "Error Estimator identity not satisfied!");
       }
       else
       {
-        libmesh_error_msg_if(std::abs(accumulated_QoI_spatially_integrated_error[0] - (-0.44809323880671958)) >= 2.e-4,
+        libmesh_error_msg_if(std::abs(time_integrated_QoI_error[0] - (-0.44809323880671958)) >= 2.e-4,
                              "Error Estimator identity not satisfied!");
 
-        libmesh_error_msg_if(std::abs(accumulated_QoI_spatially_integrated_error[1] - (-0.23895256319278346)) >= 2.e-4,
+        libmesh_error_msg_if(std::abs(time_integrated_QoI_error[1] - (-0.23895256319278346)) >= 2.e-4,
                              "Error Estimator identity not satisfied!");
       }
 #ifdef NDEBUG
diff --git a/include/solvers/adaptive_time_solver.h b/include/solvers/adaptive_time_solver.h
index f8c06f7a335..7a5eb410257 100644
--- a/include/solvers/adaptive_time_solver.h
+++ b/include/solvers/adaptive_time_solver.h
@@ -88,8 +88,11 @@ class AdaptiveTimeSolver : public FirstOrderUnsteadySolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override = 0;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override = 0;
+
   virtual Real last_completed_timestep_size() override { return completed_timestep_size; };
 
+  virtual Real get_last_deltat() override {return core_time_solver->get_last_deltat(); };
   /**
    * This method is passed on to the core_time_solver
    */
diff --git a/include/solvers/euler2_solver.h b/include/solvers/euler2_solver.h
index e73e02fce61..6de7c678aa4 100644
--- a/include/solvers/euler2_solver.h
+++ b/include/solvers/euler2_solver.h
@@ -87,6 +87,8 @@ class Euler2Solver : public FirstOrderUnsteadySolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override;
+
   /**
    * This method uses the DifferentiablePhysics'
    * element_time_derivative() and element_constraint()
diff --git a/include/solvers/euler_solver.h b/include/solvers/euler_solver.h
index c1538d9f243..da42a342dbd 100644
--- a/include/solvers/euler_solver.h
+++ b/include/solvers/euler_solver.h
@@ -109,6 +109,8 @@ class EulerSolver : public FirstOrderUnsteadySolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override;
+
   /**
    * The value for the theta method to employ: 1.0 corresponds
    * to backwards Euler, 0.0 corresponds to forwards Euler,
diff --git a/include/solvers/first_order_unsteady_solver.h b/include/solvers/first_order_unsteady_solver.h
index 85367dd649b..d0274208581 100644
--- a/include/solvers/first_order_unsteady_solver.h
+++ b/include/solvers/first_order_unsteady_solver.h
@@ -96,6 +96,8 @@ class FirstOrderUnsteadySolver : public UnsteadySolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override = 0;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override = 0;
+
 protected:
 
   /**
diff --git a/include/solvers/newmark_solver.h b/include/solvers/newmark_solver.h
index 9addccd4ba9..8a676aa0f12 100644
--- a/include/solvers/newmark_solver.h
+++ b/include/solvers/newmark_solver.h
@@ -78,6 +78,10 @@ class NewmarkSolver : public SecondOrderUnsteadySolver
    */
   virtual void adjoint_advance_timestep () override;
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override;
+
+  virtual Real get_last_deltat() override {libmesh_not_implemented();};
+
   /**
    * This method uses the specified initial displacement and velocity
    * to compute the initial acceleration \f$a_0\f$.
diff --git a/include/solvers/second_order_unsteady_solver.h b/include/solvers/second_order_unsteady_solver.h
index 56b8dfd5731..3b937220a0d 100644
--- a/include/solvers/second_order_unsteady_solver.h
+++ b/include/solvers/second_order_unsteady_solver.h
@@ -85,6 +85,10 @@ class SecondOrderUnsteadySolver : public UnsteadySolver
    */
   virtual void retrieve_timestep () override;
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override = 0;
+
+  virtual Real get_last_deltat() = 0;
+
   /**
    * Specify non-zero initial velocity. Should be called before solve().
    * The function value f and its gradient g are user-provided cloneable functors.
diff --git a/include/solvers/steady_solver.h b/include/solvers/steady_solver.h
index 85613fb8a62..c9a53ed908b 100644
--- a/include/solvers/steady_solver.h
+++ b/include/solvers/steady_solver.h
@@ -133,6 +133,10 @@ class SteadySolver : public TimeSolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override;
+
+  virtual Real get_last_deltat() override {libmesh_error_msg("Error: A steady solver is asking for a timestep size.");}
+
 protected:
 
   /**
diff --git a/include/solvers/time_solver.h b/include/solvers/time_solver.h
index 8d70c455456..2b8a4f45566 100644
--- a/include/solvers/time_solver.h
+++ b/include/solvers/time_solver.h
@@ -34,6 +34,7 @@ namespace libMesh
 class DiffContext;
 class DiffSolver;
 class DifferentiablePhysics;
+class System;
 class DifferentiableSystem;
 class ParameterVector;
 class SensitivityData;
@@ -156,6 +157,16 @@ class TimeSolver : public ReferenceCountedObject<TimeSolver>
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error);
 #endif // LIBMESH_ENABLE_AMR
 
+
+  /**
+   * A method to evaluate generic integrals of the form:
+   * int_{0}^{T} f(u,z;p) dt ~ (sum_{i=1}^{N_outer_timesteps} sum_{k=1}^{K} ( sum_{j=1}^{N_inner_steps} w_ij f_k(u^j,z^j; p^j) )
+   * The library will supply the weights w_ij and current time t_ij to the user supplied function integration_operations.
+   * Also provided will be the primal and adjoint solution(s) values, u^j, u^j-1 and z^j, z^j+1
+   * The user defined integration_operations function will evaluate w_ij * f_k(u^j,z^j; p^j).
+   */
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) = 0;
+
   /**
    * This method uses the DifferentiablePhysics
    * element_time_derivative(), element_constraint(), and
@@ -288,6 +299,12 @@ class TimeSolver : public ReferenceCountedObject<TimeSolver>
    */
   virtual Real last_completed_timestep_size();
 
+  /**
+   * Return the size of the last timestep taken to get to the current time.
+   *
+   */
+  virtual Real get_last_deltat() = 0;
+
 protected:
 
   /**
diff --git a/include/solvers/twostep_time_solver.h b/include/solvers/twostep_time_solver.h
index 8c9d55c4283..5e39f4a7f56 100644
--- a/include/solvers/twostep_time_solver.h
+++ b/include/solvers/twostep_time_solver.h
@@ -94,6 +94,7 @@ class TwostepTimeSolver : public AdaptiveTimeSolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & adjoint_refinement_error_estimator, ErrorVector & QoI_elementwise_error) override;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override;
 };
 
 
diff --git a/include/solvers/unsteady_solver.h b/include/solvers/unsteady_solver.h
index 61b420de4a9..eed7e07c53e 100644
--- a/include/solvers/unsteady_solver.h
+++ b/include/solvers/unsteady_solver.h
@@ -140,6 +140,8 @@ class UnsteadySolver : public TimeSolver
   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & /*adjoint_refinement_error_estimator*/, ErrorVector & /*QoI_elementwise_error*/) override;
 #endif // LIBMESH_ENABLE_AMR
 
+  virtual void advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations) override = 0;
+
   /**
    * This method should return the expected convergence order of the
    * (non-local) error of the time discretization scheme - e.g. 2 for the
@@ -205,6 +207,8 @@ class UnsteadySolver : public TimeSolver
     first_solve = first_solve_setting;
   }
 
+  virtual Real get_last_deltat() override {return last_deltat;};
+
 protected:
 
   /**
diff --git a/src/solvers/euler2_solver.C b/src/solvers/euler2_solver.C
index 424685c2a80..e2929782b74 100644
--- a/src/solvers/euler2_solver.C
+++ b/src/solvers/euler2_solver.C
@@ -482,4 +482,35 @@ void Euler2Solver::integrate_adjoint_refinement_error_estimate(AdjointRefinement
 }
 #endif // LIBMESH_ENABLE_AMR
 
+void Euler2Solver::advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations)
+{
+  // Enables the user to evaluate (1 - theta)*f(u_i-1) + theta*f(u_i)
+
+  // For correctness and adjoint consistency reasons we currently only support Backward Euler
+  if(theta != 0)
+   libmesh_not_implemented();
+
+  // The mesh and solutions (primal, adjoint) are already read in.
+  // So we are ready to call the user's integration operations.
+  for (auto integration_operations_iterator:integration_operations)
+    integration_operations_iterator((1.0 - theta)*_system.deltat, dynamic_cast<System&>(_system));
+
+  Real time_left = _system.time;
+
+  // Advance to t_j+1
+  _system.time = _system.time + _system.deltat;
+
+  Real time_right = _system.time;
+
+  // Retrieve the state and adjoint vectors for the next time instant
+  retrieve_timestep();
+
+  // Mesh and solution read in. Ready to call the user's integration operations.
+  for (auto integration_operations_iterator:integration_operations)
+   integration_operations_iterator(theta*(time_right - time_left), dynamic_cast<System&>(_system));
+
+  return;
+
+}
+
 } // namespace libMesh
diff --git a/src/solvers/euler_solver.C b/src/solvers/euler_solver.C
index 9e1c6b8cf0b..620e44eb1c1 100644
--- a/src/solvers/euler_solver.C
+++ b/src/solvers/euler_solver.C
@@ -284,7 +284,7 @@ void EulerSolver::integrate_adjoint_refinement_error_estimate(AdjointRefinementE
     next_step_deltat = _system.deltat;
     _system.deltat = last_step_deltat;
 
-    // The adjoint error estimate expression for a backwards facing step
+    // The adjoint error estimate expression for a backward Euler
     // scheme needs the adjoint for the last time instant, so save the current adjoint for future use
     for (auto j : make_range(_system.n_qois()))
     {
@@ -415,4 +415,51 @@ void EulerSolver::integrate_adjoint_refinement_error_estimate(AdjointRefinementE
 }
 #endif // LIBMESH_ENABLE_AMR
 
+void EulerSolver::advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations)
+{
+  // // For EulerSolver: \int_{t_i}^{t_(i+1)} f(u(t)) dt \approx f(theta u_i + (1-theta)u_i+1) (t_i+1 - t_i)
+
+  // Currently, we only support this functionality when Backward-Euler time integration is used.
+  if (theta != 1.0)
+  libmesh_not_implemented();
+
+  // u_i will be provided by the old nonlinear solution after we advance to the next time instant.
+
+  Real time_left = _system.time;
+
+  // Advance to t_j+1
+  _system.time = _system.time + _system.deltat;
+
+  Real time_right = _system.time;
+
+  // Retrieve the state and adjoint vectors for the next time instant
+  retrieve_timestep();
+
+  // Create a new weighted solution vector (1 - theta)*u_i-1 + theta*u_i
+  std::unique_ptr<NumericVector<Number>> weighted_vector = NumericVector<Number>::build(_system.comm());
+  weighted_vector->init(_system.get_vector("_old_nonlinear_solution"));
+
+  weighted_vector->add((1.0 - theta), _system.get_vector("_old_nonlinear_solution"));
+
+  weighted_vector->add(theta, *(_system.solution));
+
+  _system.solution->swap(*weighted_vector);
+
+  std::cout<<"Weighted solution norm: "<<_system.calculate_norm(*_system.solution, 0, H1)<<std::endl;
+
+  _system.update();
+
+  // Mesh and solution read in. Ready to call the user's integration operations.
+  // Since we have already weighted the solution and old solution, we just pass 1.0
+  // as the time quadrature weight.
+  for (auto integration_operations_iterator:integration_operations)
+   integration_operations_iterator(1.0*(time_right - time_left), dynamic_cast<System&>(_system));
+
+  // Swap back
+  _system.solution->swap(*weighted_vector);
+
+  return;
+}
+
+
 } // namespace libMesh
diff --git a/src/solvers/newmark_solver.C b/src/solvers/newmark_solver.C
index 98cf34ed251..d4ad9553cd7 100644
--- a/src/solvers/newmark_solver.C
+++ b/src/solvers/newmark_solver.C
@@ -104,6 +104,11 @@ void NewmarkSolver::adjoint_advance_timestep ()
   libmesh_not_implemented();
 }
 
+void NewmarkSolver::advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> /* integration_operations */)
+{
+  libmesh_not_implemented();
+}
+
 void NewmarkSolver::compute_initial_accel()
 {
   // We need to compute the initial acceleration based off of
diff --git a/src/solvers/steady_solver.C b/src/solvers/steady_solver.C
index f15f9306d56..407e913292f 100644
--- a/src/solvers/steady_solver.C
+++ b/src/solvers/steady_solver.C
@@ -133,4 +133,17 @@ void SteadySolver::integrate_adjoint_refinement_error_estimate
 }
 #endif // LIBMESH_ENABLE_AMR
 
+void SteadySolver::advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations)
+{
+  // For a steady state solve, time quadrature weight is simply 1.0
+  Real time_quadrature_weight = 1.0;
+  for(auto integration_operations_iterator:integration_operations)
+  {
+    auto f = integration_operations_iterator;
+    f(time_quadrature_weight, dynamic_cast<System&>(_system));
+  }
+
+return;
+}
+
 } // namespace libMesh
diff --git a/src/solvers/twostep_time_solver.C b/src/solvers/twostep_time_solver.C
index f4dfb45ced8..a5e56400b71 100644
--- a/src/solvers/twostep_time_solver.C
+++ b/src/solvers/twostep_time_solver.C
@@ -438,4 +438,14 @@ void TwostepTimeSolver::integrate_adjoint_refinement_error_estimate(AdjointRefin
 }
 #endif // LIBMESH_ENABLE_AMR
 
+void TwostepTimeSolver::advance_postprocessing_timestep(std::vector<std::function<void(Real, System &)>> integration_operations)
+{
+  // First half timestep
+  core_time_solver->advance_postprocessing_timestep(integration_operations);
+
+  // Second half timestep
+  core_time_solver->advance_postprocessing_timestep(integration_operations);
+}
+
+
 } // namespace libMesh
diff --git a/src/solvers/unsteady_solver.C b/src/solvers/unsteady_solver.C
index f982f0c685e..d63996b389c 100644
--- a/src/solvers/unsteady_solver.C
+++ b/src/solvers/unsteady_solver.C
@@ -237,6 +237,8 @@ void UnsteadySolver::adjoint_advance_timestep ()
 
 void UnsteadySolver::retrieve_timestep()
 {
+  last_deltat = _system.deltat;
+
   // Retrieve all the stored vectors at the current time
   solution_history->retrieve(false, _system.time);