From 112b78cbd187e3896977890fed0a69b752564194 Mon Sep 17 00:00:00 2001
From: vikramvgarg <vikram.v.garg@gmail.com>
Date: Fri, 22 Apr 2022 19:00:03 -0700
Subject: [PATCH] Added overloaded advance_postprocessing_timestep method to
 TimeSolver and its 1st level derived classes.

1) The new advance_postprocessing_timestep takes a vector of std::function as arguments and provides
these functions with the right time integration quadrature.
2) advance_postprocessing_timestep manages the time incrementation, retrieval of solutions and meshes.
---
 include/solvers/steady_solver.h   |  2 ++
 include/solvers/time_solver.h     | 11 +++++++++++
 include/solvers/unsteady_solver.h |  2 ++
 src/solvers/steady_solver.C       | 13 +++++++++++++
 4 files changed, 28 insertions(+)

diff --git a/include/solvers/steady_solver.h b/include/solvers/steady_solver.h
index 85613fb8a62..8b98b50f0c7 100644
--- a/include/solvers/steady_solver.h
+++ b/include/solvers/steady_solver.h
@@ -133,6 +133,8 @@ 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;
+
 protected:
 
   /**
diff --git a/include/solvers/time_solver.h b/include/solvers/time_solver.h
index 8d70c455456..4d9d7bf10b4 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
diff --git a/include/solvers/unsteady_solver.h b/include/solvers/unsteady_solver.h
index 61b420de4a9..93b5db76b00 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
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