Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_BackwardEuler_PseudoTransient_SA.cpp
Go to the documentation of this file.
1 //@HEADER
2 // *****************************************************************************
3 // Tempus: Time Integration and Sensitivity Analysis Package
4 //
5 // Copyright 2017 NTESS and the Tempus contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 //@HEADER
9 
12 #include "Teuchos_TimeMonitor.hpp"
13 
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp"
16 #include "Tempus_IntegratorPseudoTransientAdjointSensitivity.hpp"
17 
18 #include "Thyra_VectorStdOps.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
20 
21 #include "../TestModels/SteadyQuadraticModel.hpp"
22 
23 #include "Thyra_DefaultMultiVectorProductVector.hpp"
24 #include "Thyra_DefaultProductVector.hpp"
25 
26 #include <vector>
27 #include <fstream>
28 
29 namespace Tempus_Test {
30 
31 using Teuchos::getParametersFromXmlFile;
33 using Teuchos::RCP;
34 using Teuchos::sublist;
35 
36 // ************************************************************
37 // ************************************************************
38 void test_pseudotransient_fsa(const bool use_dfdp_as_tangent,
39  Teuchos::FancyOStream &out, bool &success)
40 {
41  // Read params from .xml file
42  RCP<ParameterList> pList =
43  getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
44 
45  // Setup the SteadyQuadraticModel
46  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
47  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
50 
51  // Setup sensitivities
52  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
53  ParameterList &sens_pl = pl->sublist("Sensitivities");
54  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
55  sens_pl.set("Reuse State Linear Solver", true);
56  sens_pl.set("Force W Update", true); // Have to do this because for this
57  // model the solver seems to be overwriting the matrix
58 
59  // Setup the Integrator
61  Tempus::createIntegratorPseudoTransientForwardSensitivity<double>(pl,
62  model);
63 
64  // Integrate to timeMax
65  bool integratorStatus = integrator->advanceTime();
66  TEST_ASSERT(integratorStatus);
67 
68  // Test if at 'Final Time'
69  double time = integrator->getTime();
70  double timeFinal = pl->sublist("Default Integrator")
71  .sublist("Time Step Control")
72  .get<double>("Final Time");
73  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
74 
75  // Time-integrated solution and the exact solution
76  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
77  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDxDp();
78  const double x = Thyra::get_ele(*x_vec, 0);
79  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
80  const double x_exact = model->getSteadyStateSolution();
81  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
82 
83  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
84  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
85 }
86 
87 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA)
88 {
89  test_pseudotransient_fsa(false, out, success);
90 }
91 
92 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_Tangent)
93 {
94  test_pseudotransient_fsa(true, out, success);
95 }
96 
97 // ************************************************************
98 // ************************************************************
99 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_ASA)
100 {
101  // Read params from .xml file
102  RCP<ParameterList> pList =
103  getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
104 
105  // Setup the SteadyQuadraticModel
106  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
109 
110  // Setup sensitivities
111  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
112  // ParameterList& sens_pl = pl->sublist("Sensitivities");
113 
114  // Setup the Integrator
116  Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl, model);
117 
118  // Integrate to timeMax
119  bool integratorStatus = integrator->advanceTime();
120  TEST_ASSERT(integratorStatus);
121 
122  // Test if at 'Final Time'
123  double time = integrator->getTime();
124  double timeFinal = pl->sublist("Default Integrator")
125  .sublist("Time Step Control")
126  .get<double>("Final Time");
127  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
128 
129  // Time-integrated solution and the exact solution using the fact that g = x
130  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
131  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
132  const double x = Thyra::get_ele(*x_vec, 0);
133  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
134  const double x_exact = model->getSteadyStateSolution();
135  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
136 
137  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
138  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
139 }
140 
141 } // namespace Tempus_Test
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
#define TEST_ASSERT(v1)
T * get() const
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Simple quadratic equation with a stable steady-state.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")