Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
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_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);
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
84TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA)
85{
86 test_pseudotransient_fsa(false, out, success);
87}
88
89TEUCHOS_UNIT_TEST(BackwardEuler, SteadyQuadratic_PseudoTransient_FSA_Tangent)
90{
91 test_pseudotransient_fsa(true, out, success);
92}
93
94// ************************************************************
95// ************************************************************
96TEUCHOS_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);
104 RCP<SteadyQuadraticModel<double> > model =
105 Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
106
107 // Setup sensitivities
108 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
109 //ParameterList& sens_pl = pl->sublist("Sensitivities");
110
111 // Setup the Integrator
112 RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<double> > 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
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)