From 201bed24c016b6bfccf31d43891be369510b2c88 Mon Sep 17 00:00:00 2001 From: "Irina K. Tezaur" Date: Tue, 2 Jun 2020 14:21:06 -0700 Subject: [PATCH] Piro: bug fix required for transient forward sensitivities to work in Albany and extension of observer to observe dx/dp (#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 --- _config.yml | 1 + packages/piro/src/Piro_ObserverBase.hpp | 84 ++++++++++++++++++- ...ToTempusIntegrationObserverAdapter_Def.hpp | 81 ++++++++++++------ .../piro/src/Piro_TransientSolver_Def.hpp | 17 ++++ ...empusSolver_SensitivitySinCosUnitTests.cpp | 29 +++++++ ...Piro_TempusSolver_SensitivityUnitTests.cpp | 53 ++++++++---- packages/piro/test/Piro_Test_MockObserver.hpp | 81 ++++++++++++++++++ 7 files changed, 303 insertions(+), 43 deletions(-) create mode 100644 _config.yml diff --git a/_config.yml b/_config.yml new file mode 100644 index 000000000000..18854876c67f --- /dev/null +++ b/_config.yml @@ -0,0 +1 @@ +theme: jekyll-theme-midnight \ No newline at end of file diff --git a/packages/piro/src/Piro_ObserverBase.hpp b/packages/piro/src/Piro_ObserverBase.hpp index 5e823eedcc3d..d6b98917fee1 100644 --- a/packages/piro/src/Piro_ObserverBase.hpp +++ b/packages/piro/src/Piro_ObserverBase.hpp @@ -46,6 +46,7 @@ #include #include "Thyra_VectorBase.hpp" +#include "Thyra_DefaultMultiVectorProductVector.hpp" namespace Piro { @@ -55,25 +56,51 @@ class ObserverBase { virtual void observeSolution( const Thyra::VectorBase &solution); + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp); + virtual void observeSolution( const Thyra::VectorBase &solution, const Scalar stamp); + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Scalar stamp); + virtual void observeSolution( const Thyra::VectorBase &solution, const Thyra::VectorBase &solution_dot, const Scalar stamp); + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Thyra::VectorBase &solution_dot, + const Scalar stamp); + virtual void observeSolution( const Thyra::VectorBase &solution, const Thyra::VectorBase &solution_dot, const Thyra::VectorBase &solution_dotdot, const Scalar stamp); - + + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Thyra::VectorBase &solution_dot, + const Thyra::VectorBase &solution_dotdot, + const Scalar stamp); virtual void observeSolution( const Thyra::MultiVectorBase &solution, Scalar time); + virtual void observeSolution( + const Thyra::MultiVectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + Scalar time); + virtual void parameterChanged( const std::string& param); @@ -92,6 +119,25 @@ template void ObserverBase::observeSolution( const Thyra::VectorBase &/*solution*/, + const Thyra::MultiVectorBase &/*solution_dxdp*/) +{ + // Nothing to do by default +} + +template +void +ObserverBase::observeSolution( + const Thyra::VectorBase &/*solution*/, + const Scalar /*stamp*/) +{ + // Nothing to do by default +} + +template +void +ObserverBase::observeSolution( + const Thyra::VectorBase &/*solution*/, + const Thyra::MultiVectorBase &/*solution_dxdp*/, const Scalar /*stamp*/) { // Nothing to do by default @@ -107,6 +153,17 @@ ObserverBase::observeSolution( // Nothing to do by default } +template +void +ObserverBase::observeSolution( + const Thyra::VectorBase &/*solution*/, + const Thyra::MultiVectorBase &/*solution_dxdp*/, + const Thyra::VectorBase &/*solution_dot*/, + const Scalar /*stamp*/) +{ + // Nothing to do by default +} + template void ObserverBase::observeSolution( @@ -121,7 +178,30 @@ ObserverBase::observeSolution( template void ObserverBase::observeSolution( - const Thyra::MultiVectorBase &solution, Scalar time) + const Thyra::VectorBase &/*solution*/, + const Thyra::MultiVectorBase &/*solution_dxdp*/, + const Thyra::VectorBase &/*solution_dot*/, + const Thyra::VectorBase &/*solution_dotdot*/, + const Scalar /*stamp*/) +{ + // Nothing to do by default +} + +template +void +ObserverBase::observeSolution( + const Thyra::MultiVectorBase &/*solution*/, + Scalar /*time*/) +{ + // Nothing to do by default +} + +template +void +ObserverBase::observeSolution( + const Thyra::MultiVectorBase &/*solution*/, + const Thyra::MultiVectorBase &/*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 index 730832da5b17..f38ef7e53990 100644 --- 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 @@ -225,25 +226,33 @@ Piro::ObserverToTempusIntegrationObserverAdapter::observeTimeStep() //Get solution typedef Thyra::DefaultMultiVectorProductVector DMVPV; - Teuchos::RCP> x = solutionHistory_->getCurrentState()->getX(); - Teuchos::RCP X = Teuchos::rcp_dynamic_cast(x); - //IKT, 5/12/2020: getX() returns a vector containing the sensitivities for the case of sensitivity calculations for sensitivity integrator. + Teuchos::RCP> x = solutionHistory_->getCurrentState()->getX(); + Teuchos::RCP X = Teuchos::rcp_dynamic_cast(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> 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> solution_dxdp = X->getMultiVector()->subView(rng); - //} - + Teuchos::RCP> solution = (sens_method_ == NONE) ? x : X->getNonconstMultiVector()->col(0); + Teuchos::RCP > 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> solution_dxdp = solution_dxdp_mv->col(np); + Teuchos::Range1D range; + RTOpPack::ConstSubVectorView 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> xdot = solutionHistory_->getCurrentState()->getXDot(); Teuchos::RCP XDot = Teuchos::rcp_dynamic_cast(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> solution_dot = (sens_method_ == NONE) ? xdot : XDot->getMultiVector()->col(0); const Scalar scalar_time = solutionHistory_->getCurrentState()->getTime(); typedef typename Teuchos::ScalarTraits::magnitudeType StampScalar; @@ -252,21 +261,41 @@ Piro::ObserverToTempusIntegrationObserverAdapter::observeTimeStep() //Get solution_dotdot Teuchos::RCP> xdotdot = solutionHistory_->getCurrentState()->getXDotDot(); Teuchos::RCP XDotDot = Teuchos::rcp_dynamic_cast(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> 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> 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> 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 index 79cf1e45247f..18c393135d1f 100644 --- a/packages/piro/src/Piro_TransientSolver_Def.hpp +++ b/packages/piro/src/Piro_TransientSolver_Def.hpp @@ -58,6 +58,8 @@ #include #include +//#define DEBUG_OUTPUT + template Piro::TransientSolver::TransientSolver( const Teuchos::RCP > &model) : @@ -402,6 +404,21 @@ Piro::TransientSolver::evalConvergedModelResponsesAndSensitivities( { //Get dxdp_mv from Tempus::ForwardIntegratorSensitivity class const RCP > 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> dxdp = dxdp_mv->col(i); +#ifdef DEBUG_OUTPUT + *out_ << "\n*** Piro::TransientSolver dxdp for p = " << i << " ***\n"; + Teuchos::Range1D range; + RTOpPack::ConstSubVectorView 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 index feeb5d28b0e5..6e986fde4c82 100644 --- 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 > observer(new MockObserver); std::vector StepSize; std::vector 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 > solutionHistory = integrator->getSolutionHistory(); + const RCP > timeStepControl = integrator->getTimeStepControl(); + const Teuchos::RCP > tempusObserver + = Teuchos::rcp(new ObserverToTempusIntegrationObserverAdapter(solutionHistory, timeStepControl, observer, false, false, sens_method)); + integrator->setObserver(tempusObserver); const RCP > stepSolver = Teuchos::null; @@ -248,6 +255,28 @@ void test_sincos_fsa(const bool use_combined_method, } ftmp.close(); } + + const RCP > solution = + observer->lastSolution(); + const RCP > 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> DxDp_vec = DxDp->col(np); + Teuchos::RCP> solution_dxdp_vec = solution_dxdp->col(np); + TEST_COMPARE_FLOATING_ARRAYS( + arrayFromVector(*solution_dxdp_vec), + arrayFromVector(*DxDp_vec), + tol); + } + // Calculate the error RCP > xdiff = x->clone_v(); RCP > DxDpdiff = DxDp->clone_mv(); diff --git a/packages/piro/test/Piro_TempusSolver_SensitivityUnitTests.cpp b/packages/piro/test/Piro_TempusSolver_SensitivityUnitTests.cpp index 2081bca3b9d1..af6e00ba11cf 100644 --- a/packages/piro/test/Piro_TempusSolver_SensitivityUnitTests.cpp +++ b/packages/piro/test/Piro_TempusSolver_SensitivityUnitTests.cpp @@ -81,6 +81,7 @@ using namespace Teuchos; using namespace Piro; using namespace Piro::Test; +//#define DEBUG_OUTPUT namespace Thyra { typedef ModelEvaluatorBase MEB; @@ -491,26 +492,48 @@ TEUCHOS_UNIT_TEST(Piro_TempusSolver, ObserveInitialConditionWhenSensitivitiesReq const Thyra::MEB::InArgs inArgs = solver->getNominalValues(); Thyra::MEB::OutArgs outArgs = solver->createOutArgs(); - { - const int solutionResponseIndex = solver->Ng() - 1; - const int parameterIndex = 0; - const Thyra::MEB::Derivative dxdp_deriv = + + const int solutionResponseIndex = solver->Ng() - 1; + const int parameterIndex = 0; + const Thyra::MEB::Derivative dxdp_deriv = Thyra::create_DgDp_mv(*solver, solutionResponseIndex, parameterIndex, Thyra::MEB::DERIV_MV_JACOBIAN_FORM); - const RCP > dxdp = dxdp_deriv.getMultiVector(); + const RCP > dxdp = dxdp_deriv.getMultiVector(); outArgs.set_DgDp(solutionResponseIndex, parameterIndex, dxdp_deriv); - } solver->evalModel(inArgs, outArgs); - { - const RCP > solution = - observer->lastSolution(); - - const RCP > initialCondition = - model->getNominalValues().get_x(); - + const RCP > solution = + observer->lastSolution(); + + const RCP > solution_dxdp = + observer->lastSolution_dxdp(); + + const RCP > initialCondition = + model->getNominalValues().get_x(); + + TEST_COMPARE_FLOATING_ARRAYS( + arrayFromVector(*solution), + arrayFromVector(*initialCondition), + tol); + + //Test observer output of lastSolution_dxdp + for (int np = 0; np < dxdp->domain()->dim(); np++) { + Teuchos::RCP> solution_dxdp_vec = solution_dxdp->col(np); + Teuchos::RCP> zero_vec = + Thyra::createMember(solution_dxdp_vec->space()); + Thyra::put_scalar(0.0, zero_vec.ptr()); +#ifdef DEBUG_OUTPUT + auto out_ =Teuchos::VerboseObjectBase::getDefaultOStream(); + *out_ << "\n*** Piro::TransientSolver dxdp for p = " << np << " ***\n"; + Teuchos::Range1D range; + RTOpPack::ConstSubVectorView dxdpv; + solution_dxdp_vec->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 = " << np << " ***\n"; +#endif TEST_COMPARE_FLOATING_ARRAYS( - arrayFromVector(*solution), - arrayFromVector(*initialCondition), + arrayFromVector(*solution_dxdp->col(np)), + arrayFromVector(*zero_vec), tol); } diff --git a/packages/piro/test/Piro_Test_MockObserver.hpp b/packages/piro/test/Piro_Test_MockObserver.hpp index 2601d9dca799..4e05e14368c4 100644 --- a/packages/piro/test/Piro_Test_MockObserver.hpp +++ b/packages/piro/test/Piro_Test_MockObserver.hpp @@ -61,25 +61,51 @@ class MockObserver : public ObserverBase { virtual void observeSolution( const Thyra::VectorBase &solution); + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp); + virtual void observeSolution( const Thyra::VectorBase &solution, const Scalar stamp); + + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Scalar stamp); virtual void observeSolution( const Thyra::VectorBase &solution, const Thyra::VectorBase &solution_dot, const Scalar stamp); + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Thyra::VectorBase &solution_dot, + const Scalar stamp); + virtual void observeSolution( const Thyra::VectorBase &solution, const Thyra::VectorBase &solution_dot, const Thyra::VectorBase &solution_dotdot, const Scalar stamp); + + virtual void observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Thyra::VectorBase &solution_dot, + const Thyra::VectorBase &solution_dotdot, + const Scalar stamp); Teuchos::RCP > lastSolution() const { return lastSolution_; } + Teuchos::RCP > lastSolution_dxdp() const { + return lastSolution_dxdp_; + } + Teuchos::RCP > lastSolution_dot() const { return lastSolution_dot_; } @@ -94,6 +120,7 @@ class MockObserver : public ObserverBase { private: Teuchos::RCP > lastSolution_; + Teuchos::RCP > lastSolution_dxdp_; Teuchos::RCP > lastSolution_dot_; Teuchos::RCP > lastSolution_dotdot_; Scalar lastStamp_; @@ -103,6 +130,7 @@ class MockObserver : public ObserverBase { template MockObserver::MockObserver() : lastSolution_(Teuchos::null), + lastSolution_dxdp_(Teuchos::null), lastSolution_dot_(Teuchos::null), lastSolution_dotdot_(Teuchos::null), lastStamp_(Scalar()) @@ -116,6 +144,17 @@ MockObserver::observeSolution( lastSolution_ = solution.clone_v(); } +template +void +MockObserver::observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp) +{ + lastSolution_ = solution.clone_v(); + lastSolution_dxdp_ = solution_dxdp.clone_mv(); +} + + template void MockObserver::observeSolution( @@ -126,15 +165,55 @@ MockObserver::observeSolution( lastStamp_ = stamp; } +template +void +MockObserver::observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Scalar stamp) +{ + lastSolution_ = solution.clone_v(); + lastSolution_dxdp_ = solution_dxdp.clone_mv(); + lastStamp_ = stamp; +} + +template +void +MockObserver::observeSolution( + const Thyra::VectorBase &solution, + const Thyra::VectorBase &solution_dot, + const Scalar stamp) +{ + lastSolution_ = solution.clone_v(); + lastSolution_dot_ = solution_dot.clone_v(); + lastStamp_ = stamp; +} + +template +void +MockObserver::observeSolution( + const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, + const Thyra::VectorBase &solution_dot, + const Scalar stamp) +{ + lastSolution_ = solution.clone_v(); + lastSolution_dxdp_ = solution_dxdp.clone_mv(); + lastSolution_dot_ = solution_dot.clone_v(); + lastStamp_ = stamp; +} + template void MockObserver::observeSolution( const Thyra::VectorBase &solution, const Thyra::VectorBase &solution_dot, + const Thyra::VectorBase &solution_dotdot, const Scalar stamp) { lastSolution_ = solution.clone_v(); lastSolution_dot_ = solution_dot.clone_v(); + lastSolution_dotdot_ = solution_dotdot.clone_v(); lastStamp_ = stamp; } @@ -142,11 +221,13 @@ template void MockObserver::observeSolution( const Thyra::VectorBase &solution, + const Thyra::MultiVectorBase &solution_dxdp, const Thyra::VectorBase &solution_dot, const Thyra::VectorBase &solution_dotdot, const Scalar stamp) { lastSolution_ = solution.clone_v(); + lastSolution_dxdp_ = solution_dxdp.clone_mv(); lastSolution_dot_ = solution_dot.clone_v(); lastSolution_dotdot_ = solution_dotdot.clone_v(); lastStamp_ = stamp;