Piro: bug fix required for transient forward sensitivities to work in…

… Albany and extension of observer to observe dx/dp (trilinos#7456) * Piro: bug fix in transient sensitivities. * Piro: bug fix to transient sensitivities. * Piro: creating observeSolution routines for observing ProductVector computing solution and sensitivities. * Piro: adding some optional debug output useful in verifying transient forward sensitivities. * Piro: fixing syntax error. * Piro: changing design for observing dx/dp. * Piro: adding some additional tests for the Tempus observer in Piro. * Piro: removing debug output. * Piro: adding some optional debug output. * Piro: adding more tests for new observer routines involving dx/dp. Co-authored-by: trilinos-autotester <[email protected]>
tcad-charon · Jun 2, 2020 · 201bed2 · 201bed2
1 parent 7ed2016
commit 201bed2
Show file tree

Hide file tree

Showing 7 changed files with 303 additions and 43 deletions.
diff --git a/_config.yml b/_config.yml
@@ -0,0 +1 @@
+theme: jekyll-theme-midnight
diff --git a/packages/piro/src/Piro_ObserverBase.hpp b/packages/piro/src/Piro_ObserverBase.hpp
@@ -46,6 +46,7 @@
 #include <string>
 
 #include "Thyra_VectorBase.hpp"
+#include "Thyra_DefaultMultiVectorProductVector.hpp"
 
 namespace Piro {
 
@@ -55,25 +56,51 @@ class ObserverBase {
   virtual void observeSolution(
       const Thyra::VectorBase<Scalar> &solution);
 
+  virtual void observeSolution(
+      const Thyra::VectorBase<Scalar> &solution,
+      const Thyra::MultiVectorBase<Scalar> &solution_dxdp); 
+
   virtual void observeSolution(
       const Thyra::VectorBase<Scalar> &solution,
       const Scalar stamp);
 
+  virtual void observeSolution(
+      const Thyra::VectorBase<Scalar> &solution,
+      const Thyra::MultiVectorBase<Scalar> &solution_dxdp,
+      const Scalar stamp);
+
   virtual void observeSolution(
       const Thyra::VectorBase<Scalar> &solution,
       const Thyra::VectorBase<Scalar> &solution_dot,
       const Scalar stamp);
 
+  virtual void observeSolution(
+      const Thyra::VectorBase<Scalar> &solution,
+      const Thyra::MultiVectorBase<Scalar> &solution_dxdp,
+      const Thyra::VectorBase<Scalar> &solution_dot,
+      const Scalar stamp);
+
   virtual void observeSolution(
       const Thyra::VectorBase<Scalar> &solution,
       const Thyra::VectorBase<Scalar> &solution_dot,
       const Thyra::VectorBase<Scalar> &solution_dotdot,
       const Scalar stamp);
-
+
+  virtual void observeSolution(
+      const Thyra::VectorBase<Scalar> &solution,
+      const Thyra::MultiVectorBase<Scalar> &solution_dxdp,
+      const Thyra::VectorBase<Scalar> &solution_dot,
+      const Thyra::VectorBase<Scalar> &solution_dotdot,
+      const Scalar stamp);
 
   virtual void observeSolution(
       const Thyra::MultiVectorBase<Scalar> &solution, Scalar time);
 
+  virtual void observeSolution(
+      const Thyra::MultiVectorBase<Scalar> &solution, 
+      const Thyra::MultiVectorBase<Scalar> &solution_dxdp,
+      Scalar time);
+
   virtual void parameterChanged(
       const std::string& param);
 
@@ -92,6 +119,25 @@ template <typename Scalar>
 void
 ObserverBase<Scalar>::observeSolution(
     const Thyra::VectorBase<Scalar> &/*solution*/,
+    const Thyra::MultiVectorBase<Scalar> &/*solution_dxdp*/)
+{
+  // Nothing to do by default
+}
+
+template <typename Scalar>
+void
+ObserverBase<Scalar>::observeSolution(
+    const Thyra::VectorBase<Scalar> &/*solution*/,
+    const Scalar /*stamp*/)
+{
+  // Nothing to do by default
+}
+
+template <typename Scalar>
+void
+ObserverBase<Scalar>::observeSolution(
+    const Thyra::VectorBase<Scalar> &/*solution*/,
+    const Thyra::MultiVectorBase<Scalar> &/*solution_dxdp*/,
     const Scalar /*stamp*/)
 {
   // Nothing to do by default
@@ -107,6 +153,17 @@ ObserverBase<Scalar>::observeSolution(
   // Nothing to do by default
 }
 
+template <typename Scalar>
+void
+ObserverBase<Scalar>::observeSolution(
+    const Thyra::VectorBase<Scalar> &/*solution*/,
+    const Thyra::MultiVectorBase<Scalar> &/*solution_dxdp*/,
+    const Thyra::VectorBase<Scalar> &/*solution_dot*/,
+    const Scalar /*stamp*/)
+{
+  // Nothing to do by default
+}
+
 template <typename Scalar>
 void
 ObserverBase<Scalar>::observeSolution(
@@ -121,7 +178,30 @@ ObserverBase<Scalar>::observeSolution(
 template <typename Scalar>
 void
 ObserverBase<Scalar>::observeSolution(
-      const Thyra::MultiVectorBase<Scalar> &solution, Scalar time)
+    const Thyra::VectorBase<Scalar> &/*solution*/,
+    const Thyra::MultiVectorBase<Scalar> &/*solution_dxdp*/,
+    const Thyra::VectorBase<Scalar> &/*solution_dot*/,
+    const Thyra::VectorBase<Scalar> &/*solution_dotdot*/,
+    const Scalar /*stamp*/)
+{
+  // Nothing to do by default
+}
+
+template <typename Scalar>
+void
+ObserverBase<Scalar>::observeSolution(
+      const Thyra::MultiVectorBase<Scalar> &/*solution*/, 
+      Scalar /*time*/)
+{
+  // Nothing to do by default
+}
+
+template <typename Scalar>
+void
+ObserverBase<Scalar>::observeSolution(
+      const Thyra::MultiVectorBase<Scalar> &/*solution*/, 
+      const Thyra::MultiVectorBase<Scalar> &/*solution_dxdp*/,
+      Scalar /*time*/)
 {
   // Nothing to do by default
 }

diff --git a/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp b/packages/piro/src/Piro_ObserverToTempusIntegrationObserverAdapter_Def.hpp
@@ -44,6 +44,7 @@
 
 #include "Teuchos_Ptr.hpp"
 
+//#define DEBUG_OUTPUT
 
 // Constructor
 template <typename Scalar>
@@ -225,25 +226,33 @@ Piro::ObserverToTempusIntegrationObserverAdapter<Scalar>::observeTimeStep()
 
   //Get solution
   typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
-  Teuchos::RCP<const Thyra::VectorBase<Scalar>> x = solutionHistory_->getCurrentState()->getX(); 
-  Teuchos::RCP<const DMVPV> X = Teuchos::rcp_dynamic_cast<const DMVPV>(x);
-  //IKT, 5/12/2020: getX() returns a vector containing the sensitivities for the case of sensitivity calculations for sensitivity integrator. 
+  Teuchos::RCP<Thyra::VectorBase<Scalar>> x = solutionHistory_->getCurrentState()->getX(); 
+  Teuchos::RCP<DMVPV> X = Teuchos::rcp_dynamic_cast<DMVPV>(x);
+  //IKT, 5/12/2020: getX() returns a vector containing the sensitivities for the 
+  //case of sensitivity calculations for sensitivity integrator. 
   //Hence, we extract the first column, which is the solution.
-  Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution = (sens_method_ == NONE) ? x : X->getMultiVector()->col(0);
-  //IKT FIXME? Extract the remaining columns which would represent dxdp for the forward sensitivities?
-  //This is if we want to observe them/write them to Exodus file in Albany.
-  //if (sens_method == FORWARD) {
-    //const int num_param = X->getMultiVector()->domain()->dim()-1;
-    //const Teuchos::Range1D rng(1,num_param);
-    //Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dxdp = X->getMultiVector()->subView(rng);
-  //}
-
+  Teuchos::RCP<Thyra::VectorBase<Scalar>> solution = (sens_method_ == NONE) ? x : X->getNonconstMultiVector()->col(0);
+  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > solution_dxdp_mv = Teuchos::null; 
+  if (sens_method_ == FORWARD) {
+    const int num_param = X->getMultiVector()->domain()->dim()-1;
+    const Teuchos::Range1D rng(1,num_param);
+    solution_dxdp_mv = X->getNonconstMultiVector()->subView(rng);
+#ifdef DEBUG_OUTPUT
+    for (int np = 0; np < num_param; np++) {
+      *out_ << "\n*** Piro::ObserverToTempusIntegrationObserverAdapter dxdp" << np << " ***\n";
+      Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dxdp = solution_dxdp_mv->col(np);
+      Teuchos::Range1D range;
+      RTOpPack::ConstSubVectorView<Scalar> dxdpv;
+      solution_dxdp->acquireDetachedView(range, &dxdpv);
+      auto dxdpa = dxdpv.values();
+      for (auto i = 0; i < dxdpa.size(); ++i) *out_ << dxdpa[i] << " ";
+      *out_ << "\n*** Piro::ObserverToTempusIntegrationObserverAdapter dxdp" << np << " ***\n";
+    }
+#endif
+  }
   //Get solution_dot 
   Teuchos::RCP<const Thyra::VectorBase<Scalar>> xdot = solutionHistory_->getCurrentState()->getXDot(); 
   Teuchos::RCP<const DMVPV> XDot = Teuchos::rcp_dynamic_cast<const DMVPV>(xdot);
-  //IKT, 5/11/2020: getXDot() returns a vector containing the sensitivities for the case of sensitivity calculations for sentivity integrator 
-  //Hence, we extract the first column, which is the solution_dot.
-  Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dot = (sens_method_ == NONE) ? xdot : XDot->getMultiVector()->col(0);
 
   const Scalar scalar_time = solutionHistory_->getCurrentState()->getTime();
   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType StampScalar;
@@ -252,21 +261,41 @@ Piro::ObserverToTempusIntegrationObserverAdapter<Scalar>::observeTimeStep()
   //Get solution_dotdot 
   Teuchos::RCP<const Thyra::VectorBase<Scalar>> xdotdot = solutionHistory_->getCurrentState()->getXDotDot(); 
   Teuchos::RCP<const DMVPV> XDotDot = Teuchos::rcp_dynamic_cast<const DMVPV>(xdotdot);
-  //IKT, 5/11/2020: getXDotDot() returns a vector containing the sensitivities for the case of sensitivity calculations. 
-  //Hence, we extract the first column, which is the solution_dotdot.
-  Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dotdot = (sens_method_ == NONE) ? xdot : XDotDot->getMultiVector()->col(0);
-  if (Teuchos::nonnull(solution_dot))
-  {
+  if (Teuchos::nonnull(xdot)) {
+    //IKT, 5/11/2020: getXDot() returns a vector containing the sensitivities for 
+    //the case of sensitivity calculations for sentivity integrator 
+    //Hence, we extract the first column, which is the solution_dot.
+    Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dot 
+          = (sens_method_ == NONE) ? xdot : XDot->getMultiVector()->col(0);
     if (supports_x_dotdot_) {
-
-      wrappedObserver_->observeSolution(*solution, *solution_dot, *solution_dotdot, time);
+      //IKT, 5/11/2020: getXDotDot() returns a vector containing the sensitivities for 
+      //the case of sensitivity calculations. 
+      //Hence, we extract the first column, which is the solution_dotdot.
+      Teuchos::RCP<const Thyra::VectorBase<Scalar>> solution_dotdot 
+          = (sens_method_ == NONE) ? xdot : XDotDot->getMultiVector()->col(0);
+      if (solution_dxdp_mv != Teuchos::null) {
+        wrappedObserver_->observeSolution(*solution, *solution_dxdp_mv, *solution_dot, *solution_dotdot, time);
+      }
+      else {
+        wrappedObserver_->observeSolution(*solution, *solution_dot, *solution_dotdot, time);
+      }
+    }
+    else {
+      if (solution_dxdp_mv != Teuchos::null) {
+        wrappedObserver_->observeSolution(*solution, *solution_dxdp_mv, *solution_dot, time);
+      }
+      else {
+        wrappedObserver_->observeSolution(*solution, *solution_dot, time);
+      }
     }
-   else {
-      wrappedObserver_->observeSolution(*solution, *solution_dot, time);
-   }
   }
   else {
-    wrappedObserver_->observeSolution(*solution, time);
+    if (solution_dxdp_mv != Teuchos::null) {
+      wrappedObserver_->observeSolution(*solution, *solution_dxdp_mv, time);
+    }
+    else {
+      wrappedObserver_->observeSolution(*solution, time);
+    }
   }
   previous_dt_ = current_dt; 
 }

diff --git a/packages/piro/src/Piro_TransientSolver_Def.hpp b/packages/piro/src/Piro_TransientSolver_Def.hpp
@@ -58,6 +58,8 @@
 #include <stdexcept>
 #include <iostream>
 
+//#define DEBUG_OUTPUT
+
 template <typename Scalar>
 Piro::TransientSolver<Scalar>::TransientSolver(
   const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &model) :  
@@ -402,6 +404,21 @@ Piro::TransientSolver<Scalar>::evalConvergedModelResponsesAndSensitivities(
       {
         //Get dxdp_mv from Tempus::ForwardIntegratorSensitivity class  
         const RCP<const Thyra::MultiVectorBase<Scalar> > dxdp_mv = piroTempusIntegrator_->getDxDp();
+#ifdef DEBUG_OUTPUT
+        *out_ << "\n*** Piro::TransientSolver: num_p, num vecs in dxdp = " << num_p_ << ", " << dxdp_mv->domain()->dim() << " ***\n";
+#endif
+        for (int i=0; i < dxdp_mv->domain()->dim(); ++i) { 
+          Teuchos::RCP<const Thyra::VectorBase<Scalar>> dxdp = dxdp_mv->col(i);
+#ifdef DEBUG_OUTPUT
+          *out_ << "\n*** Piro::TransientSolver dxdp for p = " << i << " ***\n";
+          Teuchos::Range1D range;
+          RTOpPack::ConstSubVectorView<Scalar> dxdpv;
+          dxdp->acquireDetachedView(range, &dxdpv);
+          auto dxdpa = dxdpv.values();
+          for (auto j = 0; j < dxdpa.size(); ++j) *out_ << dxdpa[j] << " ";
+          *out_ << "\n*** Piro::TransientSolver dxdp for p = " << i << " ***\n";
+#endif
+        }
         //IMPORTANT REMARK: we are currently NOT using DxdotDp and DxdotdotDp in transient sensitivities!  
         //The capability to use them can be added at a later point in time, if desired. 
         //IKT, 5/10/20: throw error if dxdp_mv returned by Tempus is null.  Not sure if this can happen in practice or not...

diff --git a/packages/piro/test/Piro_TempusSolver_SensitivitySinCosUnitTests.cpp b/packages/piro/test/Piro_TempusSolver_SensitivitySinCosUnitTests.cpp
@@ -53,6 +53,7 @@
 #include "Piro_TempusIntegrator.hpp"
 
 #include "SinCosModel.hpp"
+#include "Piro_Test_MockObserver.hpp"
 
 #include "Teuchos_UnitTestHarness.hpp"
 
@@ -103,6 +104,7 @@ void test_sincos_fsa(const bool use_combined_method,
     }
   }
 
+  const RCP<MockObserver<double> > observer(new MockObserver<double>);
   std::vector<double> StepSize;
   std::vector<double> ErrorNorm;
   const int nTimeStepSizes = 7;
@@ -167,6 +169,11 @@ void test_sincos_fsa(const bool use_combined_method,
     }
     integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
                                 DxDp0, Teuchos::null, Teuchos::null);
+    const RCP<const Tempus::SolutionHistory<double> > solutionHistory = integrator->getSolutionHistory();
+    const RCP<const Tempus::TimeStepControl<double> > timeStepControl = integrator->getTimeStepControl();
+    const Teuchos::RCP<Tempus::IntegratorObserver<double> > tempusObserver 
+          = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter<double>(solutionHistory, timeStepControl, observer, false, false, sens_method));
+    integrator->setObserver(tempusObserver);
 
     const RCP<Thyra::NonlinearSolverBase<double> > stepSolver = Teuchos::null;
 
@@ -248,6 +255,28 @@ void test_sincos_fsa(const bool use_combined_method,
       }
       ftmp.close();
     }
+
+    const RCP<const Thyra::VectorBase<double> > solution =
+      observer->lastSolution();
+    const RCP<const Thyra::MultiVectorBase<double> > solution_dxdp =
+      observer->lastSolution_dxdp();
+
+    //Compare solution from observer and x to verify observer routines 
+    TEST_COMPARE_FLOATING_ARRAYS(
+      arrayFromVector(*solution),
+      arrayFromVector(*x),
+      tol);
+
+    //Compare solution_dxdp from observer and DxDp to verify observer routines 
+    for (int np = 0; np < DxDp->domain()->dim(); np++) { 
+      Teuchos::RCP<const Thyra::VectorBase<double>> DxDp_vec = DxDp->col(np);
+      Teuchos::RCP<const Thyra::VectorBase<double>> solution_dxdp_vec = solution_dxdp->col(np);
+      TEST_COMPARE_FLOATING_ARRAYS(
+        arrayFromVector(*solution_dxdp_vec),
+        arrayFromVector(*DxDp_vec),
+        tol);
+    }
+
     // Calculate the error
     RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
     RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();