Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_BDF2_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::getParametersFromXmlFile;
32 using Teuchos::RCP;
33 using Teuchos::sublist;
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_BDF2_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);
47  RCP<SteadyQuadraticModel<double> > model =
48  Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
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
59  RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<double> > integrator =
60  Tempus::createIntegratorPseudoTransientForwardSensitivity<double>(pl,
61  model);
62 
63  // Integrate to timeMax
64  bool integratorStatus = integrator->advanceTime();
65  TEST_ASSERT(integratorStatus);
66 
67  // Test if at 'Final Time'
68  double time = integrator->getTime();
69  double timeFinal = pl->sublist("Default Integrator")
70  .sublist("Time Step Control")
71  .get<double>("Final Time");
72  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
73 
74  // Time-integrated solution and the exact solution
75  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
76  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDxDp();
77  const double x = Thyra::get_ele(*x_vec, 0);
78  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
79  const double x_exact = model->getSteadyStateSolution();
80  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
81 
82  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
83  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
84 }
85 
86 TEUCHOS_UNIT_TEST(BDF2, SteadyQuadratic_PseudoTransient_FSA)
87 {
88  test_pseudotransient_fsa(false, out, success);
89 }
90 
91 TEUCHOS_UNIT_TEST(BDF2, SteadyQuadratic_PseudoTransient_FSA_Tangent)
92 {
93  test_pseudotransient_fsa(true, out, success);
94 }
95 
96 // ************************************************************
97 // ************************************************************
98 TEUCHOS_UNIT_TEST(BDF2, SteadyQuadratic_PseudoTransient_ASA)
99 {
100  // Read params from .xml file
101  RCP<ParameterList> pList =
102  getParametersFromXmlFile("Tempus_BDF2_SteadyQuadratic.xml");
103 
104  // Setup the SteadyQuadraticModel
105  RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
108 
109  // Setup sensitivities
110  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
111  // ParameterList& sens_pl = pl->sublist("Sensitivities");
112 
113  // Setup the Integrator
115  Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl, model);
116 
117  // Integrate to timeMax
118  bool integratorStatus = integrator->advanceTime();
119  TEST_ASSERT(integratorStatus);
120 
121  // Test if at 'Final Time'
122  double time = integrator->getTime();
123  double timeFinal = pl->sublist("Default Integrator")
124  .sublist("Time Step Control")
125  .get<double>("Final Time");
126  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
127 
128  // Time-integrated solution and the exact solution using the fact that g = x
129  RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
130  RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
131  const double x = Thyra::get_ele(*x_vec, 0);
132  const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
133  const double x_exact = model->getSteadyStateSolution();
134  const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
135 
136  TEST_FLOATING_EQUALITY(x, x_exact, 1.0e-6);
137  TEST_FLOATING_EQUALITY(dxdb, dxdb_exact, 1.0e-6);
138 }
139 
140 } // namespace Tempus_Test
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
#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.