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: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
11 #include "Teuchos_TimeMonitor.hpp"
12 
13 #include "Tempus_config.hpp"
14 #include "Tempus_IntegratorPseudoTransientForwardSensitivity.hpp"
15 #include "Tempus_IntegratorPseudoTransientAdjointSensitivity.hpp"
16 
17 #include "Thyra_VectorStdOps.hpp"
18 #include "Thyra_MultiVectorStdOps.hpp"
19 
20 #include "../TestModels/SteadyQuadraticModel.hpp"
21 
22 #include "Thyra_DefaultMultiVectorProductVector.hpp"
23 #include "Thyra_DefaultProductVector.hpp"
24 
25 #include <vector>
26 #include <fstream>
27 
28 namespace Tempus_Test {
29 
30 using Teuchos::RCP;
32 using Teuchos::sublist;
33 using Teuchos::getParametersFromXmlFile;
34 
35 // ************************************************************
36 // ************************************************************
37 void test_pseudotransient_fsa(const bool use_dfdp_as_tangent,
38  Teuchos::FancyOStream &out, bool &success)
39 {
40  // Read params from .xml file
41  RCP<ParameterList> pList =
42  getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
43 
44  // Setup the SteadyQuadraticModel
45  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
46  scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
49 
50  // Setup sensitivities
51  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
52  ParameterList& sens_pl = pl->sublist("Sensitivities");
53  sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
54  sens_pl.set("Reuse State Linear Solver", true);
55  sens_pl.set("Force W Update", true); // Have to do this because for this
56  // model the solver seems to be overwriting the matrix
57 
58  // Setup the Integrator
60  Tempus::integratorPseudoTransientForwardSensitivity<double>(pl, model);
61 
62  // Integrate to timeMax
63  bool integratorStatus = integrator->advanceTime();
64  TEST_ASSERT(integratorStatus);
65 
66  // Test if at 'Final Time'
67  double time = integrator->getTime();
68  double timeFinal =pl->sublist("Default Integrator")
69  .sublist("Time Step Control").get<double>("Final Time");
70  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
71 
72  // Time-integrated solution and the exact solution
73  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
74  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDxDp();
75  const double x = Thyra::get_ele(*x_vec, 0);
76  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
77  const double x_exact = model->getSteadyStateSolution();
78  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
79 
80  TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 );
81  TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 );
82 }
83 
84 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA)
85 {
86  test_pseudotransient_fsa(false, out, success);
87 }
88 
89 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_Tangent)
90 {
91  test_pseudotransient_fsa(true, out, success);
92 }
93 
94 // ************************************************************
95 // ************************************************************
96 TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_ASA)
97 {
98  // Read params from .xml file
99  RCP<ParameterList> pList =
100  getParametersFromXmlFile("Tempus_BackwardEuler_SteadyQuadratic.xml");
101 
102  // Setup the SteadyQuadraticModel
103  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
106 
107  // Setup sensitivities
108  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
109  //ParameterList& sens_pl = pl->sublist("Sensitivities");
110 
111  // Setup the Integrator
113  Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl, model);
114 
115  // Integrate to timeMax
116  bool integratorStatus = integrator->advanceTime();
117  TEST_ASSERT(integratorStatus);
118 
119  // Test if at 'Final Time'
120  double time = integrator->getTime();
121  double timeFinal =pl->sublist("Default Integrator")
122  .sublist("Time Step Control").get<double>("Final Time");
123  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
124 
125  // Time-integrated solution and the exact solution using the fact that g = x
126  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
127  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
128  const double x = Thyra::get_ele(*x_vec, 0);
129  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
130  const double x_exact = model->getSteadyStateSolution();
131  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
132 
133  TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 );
134  TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 );
135 }
136 
137 } // 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. This is a simple differential equation which h...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")