Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_DIRK_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
9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
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
28namespace Tempus_Test {
29
30using Teuchos::RCP;
31using Teuchos::ParameterList;
32using Teuchos::sublist;
33using Teuchos::getParametersFromXmlFile;
34
35// ************************************************************
36// ************************************************************
37void 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_DIRK_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, 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
85TEUCHOS_UNIT_TEST(DIRK, SteadyQuadratic_PseudoTransient_FSA)
86{
87 test_pseudotransient_fsa(false, out, success);
88}
89
90TEUCHOS_UNIT_TEST(DIRK, SteadyQuadratic_PseudoTransient_FSA_Tangent)
91{
92 test_pseudotransient_fsa(true, out, success);
93}
94
95// ************************************************************
96// ************************************************************
97TEUCHOS_UNIT_TEST(DIRK, SteadyQuadratic_PseudoTransient_ASA)
98{
99 // Read params from .xml file
100 RCP<ParameterList> pList =
101 getParametersFromXmlFile("Tempus_DIRK_SteadyQuadratic.xml");
102
103 // Setup the SteadyQuadraticModel
104 RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
105 RCP<SteadyQuadraticModel<double> > model =
106 Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
107
108 // Setup sensitivities
109 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
110 //ParameterList& sens_pl = pl->sublist("Sensitivities");
111
112 // Setup the Integrator
113 RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<double> > integrator =
114 Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl, model);
115
116 // Integrate to timeMax
117 bool integratorStatus = integrator->advanceTime();
118 TEST_ASSERT(integratorStatus);
119
120 // Test if at 'Final Time'
121 double time = integrator->getTime();
122 double timeFinal =pl->sublist("Default Integrator")
123 .sublist("Time Step Control").get<double>("Final Time");
124 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
125
126 // Time-integrated solution and the exact solution using the fact that g = x
127 RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
128 RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
129 const double x = Thyra::get_ele(*x_vec, 0);
130 const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
131 const double x_exact = model->getSteadyStateSolution();
132 const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
133
134 TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 );
135 TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 );
136}
137
138
139} // namespace Tempus_Test
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)