Tempus_BackwardEuler_PseudoTransient_SA.cpp

Go to the documentation of this file.

//@HEADER
// *****************************************************************************
//          Tempus: Time Integration and Sensitivity Analysis Package
//
// Copyright 2017 NTESS and the Tempus contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
//@HEADER
 
#include "Teuchos_UnitTestHarness.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_TimeMonitor.hpp"
 
#include "Tempus_config.hpp"
#include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp"
#include "Tempus_IntegratorPseudoTransientAdjointSensitivity.hpp"
 
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
 
#include "../TestModels/SteadyQuadraticModel.hpp"
 
#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_DefaultProductVector.hpp"
 
#include <vector>
#include <fstream>
 
namespace Tempus_Test {
 
using Teuchos::getParametersFromXmlFile;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::sublist;
 
// NOX row sum scaling for the stepper of the SteadyQuadratic parameters
void setRowSumScaling(ParameterList &pl)
{
  pl.sublist("Default Stepper")
      .sublist("Default Solver")
      .sublist("NOX")
      .sublist("Thyra Group Options")
      .set("Function Scaling", "Row Sum");
}
 
// ************************************************************
// ************************************************************
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent,
                              Teuchos::FancyOStream &out, bool &success,
                              const bool use_row_sum_scaling = false)
{
  // Read params from .xml file
  RCP<ParameterList> pList =
      getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
 
  // Setup the SteadyQuadraticModel
  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
  RCP<SteadyQuadraticModel<double> > model =
      Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
 
  // Setup sensitivities
  RCP<ParameterList> pl  = sublist(pList, "Tempus", true);
  ParameterList &sens_pl = pl->sublist("Sensitivities");
  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
  // With row sum scaling, the sensitivity integrator has its own solver
  sens_pl.set("Reuse State Linear Solver", !use_row_sum_scaling);
  sens_pl.set("Force W Update", true);  // Have to do this because for this
  // model the solver seems to be overwriting the matrix
  if (use_row_sum_scaling) setRowSumScaling(*pl);
 
  // Setup the Integrator
  RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<double> > integrator =
      Tempus::createIntegratorPseudoTransientForwardSensitivity<double>(pl,
                                                                        model);
 
  // Integrate to timeMax
  bool integratorStatus = integrator->advanceTime();
  TEST_ASSERT(integratorStatus);
 
  // Test if at 'Final Time'
  double time      = integrator->getTime();
  double timeFinal = pl->sublist("Default Integrator")
                         .sublist("Time Step Control")
                         .get<double>("Final Time");
  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
 
  // Time-integrated solution and the exact solution
  RCP<const Thyra::VectorBase<double> > x_vec         = integrator->getX();
  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDxDp();
  const double x                                      = Thyra::get_ele(*x_vec, 0);
  const double dxdb                                   = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
  const double x_exact                                = model->getSteadyStateSolution();
  const double dxdb_exact                             = model->getSteadyStateSolutionSensitivity();
 
  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
}
 
TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA)
{
  test_pseudotransient_fsa(false, out, success);
}
 
TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_Tangent)
{
  test_pseudotransient_fsa(true, out, success);
}
 
// NOX row sum scaling: the state and sensitivity integrators are created from
// the same parameters, in which the NOX solver of the state stepper stores its
// row sum scaling objects (in the state residual space).  The sensitivity
// integrator must not use them.
TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_RowSumScaling)
{
  test_pseudotransient_fsa(false, out, success, true);
}
 
// ************************************************************
// ************************************************************
void test_pseudotransient_asa(Teuchos::FancyOStream &out, bool &success,
                              const bool use_row_sum_scaling = false)
{
  // Read params from .xml file
  RCP<ParameterList> pList =
      getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
 
  // Setup the SteadyQuadraticModel
  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
  RCP<SteadyQuadraticModel<double> > model =
      Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
 
  // Setup sensitivities
  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
  // ParameterList& sens_pl = pl->sublist("Sensitivities");
  if (use_row_sum_scaling) setRowSumScaling(*pl);
 
  // Setup the Integrator.  The implicit adjoint operator does not provide the
  // row sums needed by row sum scaling, so use an explicit adjoint model: the
  // model is scalar, hence it is its own adjoint.
  RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<double> > integrator =
      use_row_sum_scaling
          ? Tempus::integratorPseudoTransientAdjointSensitivity<double>(
                pl, model, model)
          : Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl,
                                                                        model);
 
  // Integrate to timeMax
  bool integratorStatus = integrator->advanceTime();
  TEST_ASSERT(integratorStatus);
 
  // Test if at 'Final Time'
  double time      = integrator->getTime();
  double timeFinal = pl->sublist("Default Integrator")
                         .sublist("Time Step Control")
                         .get<double>("Final Time");
  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
 
  // Time-integrated solution and the exact solution using the fact that g = x
  RCP<const Thyra::VectorBase<double> > x_vec         = integrator->getX();
  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
  const double x                                      = Thyra::get_ele(*x_vec, 0);
  const double dxdb                                   = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
  const double x_exact                                = model->getSteadyStateSolution();
  const double dxdb_exact                             = model->getSteadyStateSolutionSensitivity();
 
  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
}
 
TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_ASA)
{
  test_pseudotransient_asa(out, success);
}
 
// NOX row sum scaling: the state and adjoint integrators are created from the
// same parameters, in which the NOX solver of the state stepper stores its row
// sum scaling objects (in the state residual space).  The adjoint integrator
// must not use them.
TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_ASA_RowSumScaling)
{
  test_pseudotransient_asa(out, success, true);
}
 
}  // namespace Tempus_Test