Add advance_postprocessing_timestep declarations and implementations …

…for 1st order fixed timestep solvers.

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()

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,

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:
 
   /**

src/solvers/euler2_solver.C

@@ -482,4 +482,31 @@ 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, dynamic_cast<System&>(_system));
+
+  // Advance to t_j+1
+  _system.time = _system.time + _system.deltat;
+
+  // 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, dynamic_cast<System&>(_system));
+
+  return;
+
+}
+
 } // namespace libMesh

src/solvers/euler_solver.C

@@ -415,4 +415,43 @@ 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.
+
+  // Advance to t_j+1
+  _system.time = _system.time + _system.deltat;
+
+  // 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);
+
+  // 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, dynamic_cast<System&>(_system));
+
+  // Swap back
+  _system.solution->swap(*weighted_vector);
+
+  return;
+}
+
+
 } // namespace libMesh