Merge branch 'develop' of github.com:LLNL/serac into feature/bramwell…

…/strumpack

src/serac/numerics/odes.cpp

@@ -21,6 +21,8 @@ SecondOrderODE::SecondOrderODE(int n, State&& state, const EquationSolver& solve
 
 void SecondOrderODE::SetTimestepper(const serac::TimestepMethod timestepper)
 {
+  timestepper_ = timestepper;
+
   switch (timestepper) {
     case serac::TimestepMethod::Newmark:
       second_order_ode_solver_ = std::make_unique<mfem::NewmarkSolver>();
@@ -258,6 +260,8 @@ FirstOrderODE::FirstOrderODE(int n, FirstOrderODE::State&& state, const Equation
 
 void FirstOrderODE::SetTimestepper(const serac::TimestepMethod timestepper)
 {
+  timestepper_ = timestepper;
+
   switch (timestepper) {
     case serac::TimestepMethod::BackwardEuler:
       ode_solver_ = std::make_unique<mfem::BackwardEulerSolver>();

src/serac/numerics/odes.hpp

@@ -163,6 +163,13 @@ class SecondOrderODE : public mfem::SecondOrderTimeDependentOperator {
    */
   const State& GetState() { return state_; }
 
+  /**
+   * @brief Query the timestep method for the ode solver
+   *
+   * @return The timestep method used by the underlying ode solver
+   */
+  TimestepMethod GetTimestepper() { return timestepper_; }
+
 private:
   /**
    * @brief Internal implementation used for mfem::SOTDO::Mult and mfem::SOTDO::ImplicitSolve
@@ -214,6 +221,8 @@ class SecondOrderODE : public mfem::SecondOrderTimeDependentOperator {
   mutable mfem::Vector U_plus_;
   mutable mfem::Vector dU_dt_;
   mutable mfem::Vector d2U_dt2_;
+
+  serac::TimestepMethod timestepper_;
 };
 
 /**
@@ -338,6 +347,13 @@ class FirstOrderODE : public mfem::TimeDependentOperator {
     }
   }
 
+  /**
+   * @brief Query the timestep method for the ode solver
+   *
+   * @return The timestep method used by the underlying ode solver
+   */
+  TimestepMethod GetTimestepper() { return timestepper_; }
+
 private:
   /**
    * @brief Internal implementation used for mfem::TDO::Mult and mfem::TDO::ImplicitSolve\
@@ -378,6 +394,8 @@ class FirstOrderODE : public mfem::TimeDependentOperator {
   mutable mfem::Vector U_;
   mutable mfem::Vector U_plus_;
   mutable mfem::Vector dU_dt_;
+
+  TimestepMethod timestepper_;
 };
 
 }  // namespace serac::mfem_ext

src/serac/physics/base_physics.hpp

@@ -115,7 +115,7 @@ class BasePhysics {
    * @param parameter_name The name of the parameter to generate
    *
    * @note The user is responsible for managing the lifetime of this object. It is required
-   * to exist whenever advanceTimestep, solveAdjoint, or computeSensitivity is called.
+   * to exist whenever advanceTimestep, reverseAdjointTimestep, or computeTimestepSensitivity is called.
    *
    * @note The finite element space for this object is generated from the parameter
    * discretization space (e.g. L2, H1) and the computational mesh given in the physics module constructor.
@@ -190,9 +190,9 @@ class BasePhysics {
    *
    * @return The sensitivity with respect to the parameter
    *
-   * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method.
+   * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method.
    */
-  virtual FiniteElementDual& computeSensitivity(size_t /* parameter_index */)
+  virtual FiniteElementDual& computeTimestepSensitivity(size_t /* parameter_index */)
   {
     SLIC_ERROR_ROOT(axom::fmt::format("Parameter sensitivities not enabled in physics module {}", name_));
     return *parameters_[0].sensitivity;
@@ -204,14 +204,28 @@ class BasePhysics {
    *
    * @return The sensitivity with respect to the shape displacement
    *
-   * @pre `solveAdjoint` with an appropriate adjoint load must be called prior to this method.
+   * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method.
    */
-  virtual FiniteElementDual& computeShapeSensitivity()
+  virtual FiniteElementDual& computeTimestepShapeSensitivity()
   {
     SLIC_ERROR_ROOT(axom::fmt::format("Shape sensitivities not enabled in physics module {}", name_));
     return shape_displacement_sensitivity_;
   }
 
+  /**
+   * @brief Compute the implicit sensitivity of the quantity of interest with respect to the initial condition fields
+   *
+   * @return Fields states corresponding to the sensitivities with respect to the initial condition fields
+   *
+   * @pre `reverseAdjointTimestep` with an appropriate adjoint load must be called prior to this method as many times as
+   * the forward advance is called.
+   */
+  virtual const std::unordered_map<std::string, const serac::FiniteElementDual&> computeInitialConditionSensitivity()
+  {
+    SLIC_ERROR_ROOT(axom::fmt::format("Initial condition sensitivities not enabled in physics module {}", name_));
+    return {};
+  }
+
   /**
    * @brief Advance the state variables according to the chosen time integrator
    *
@@ -227,7 +241,7 @@ class BasePhysics {
    *
    * @return The computed adjoint finite element states
    */
-  virtual const std::unordered_map<std::string, const serac::FiniteElementState&> solveAdjoint(
+  virtual const std::unordered_map<std::string, const serac::FiniteElementState&> reverseAdjointTimestep(
       std::unordered_map<std::string, const serac::FiniteElementDual&> /* adjoint_loads */,
       std::unordered_map<std::string, const serac::FiniteElementState&> /* adjoint_with_essential_boundary */ = {})
   {