Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_ExplicitRK_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_ExplicitRK_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
55 // Setup the Integrator
56 RCP<Tempus::IntegratorPseudoTransientForwardSensitivity<double> > integrator =
57 Tempus::createIntegratorPseudoTransientForwardSensitivity<double>(pl, model);
58
59 // Integrate to timeMax
60 bool integratorStatus = integrator->advanceTime();
61 TEST_ASSERT(integratorStatus);
62
63 // Test if at 'Final Time'
64 double time = integrator->getTime();
65 double timeFinal = pl->sublist("Demo Integrator")
66 .sublist("Time Step Control").get<double>("Final Time");
67 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
68
69 // Time-integrated solution and the exact solution
70 RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
71 RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDxDp();
72 const double x = Thyra::get_ele(*x_vec, 0);
73 const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
74 const double x_exact = model->getSteadyStateSolution();
75 const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
76
77 TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 );
78 TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 );
79}
80
81
82TEUCHOS_UNIT_TEST(ExplicitRK, SteadyQuadratic_PseudoTransient_FSA)
83{
84 test_pseudotransient_fsa(false, out, success);
85}
86
87TEUCHOS_UNIT_TEST(ExplicitRK, SteadyQuadratic_PseudoTransient_FSA_Tangent)
88{
89 test_pseudotransient_fsa(true, out, success);
90}
91
92// ************************************************************
93// ************************************************************
94TEUCHOS_UNIT_TEST(ExplicitRK, SteadyQuadratic_PseudoTransient_ASA)
95{
96 // Read params from .xml file
97 RCP<ParameterList> pList =
98 getParametersFromXmlFile("Tempus_ExplicitRK_SteadyQuadratic.xml");
99
100 // Setup the SteadyQuadraticModel
101 RCP<ParameterList> scm_pl = sublist(pList, "SteadyQuadraticModel", true);
102 RCP<SteadyQuadraticModel<double> > model =
103 Teuchos::rcp(new SteadyQuadraticModel<double>(scm_pl));
104
105 // Setup sensitivities
106 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
107 //ParameterList& sens_pl = pl->sublist("Sensitivities");
108
109 // Setup the Integrator
110 RCP<Tempus::IntegratorPseudoTransientAdjointSensitivity<double> > integrator =
111 Tempus::integratorPseudoTransientAdjointSensitivity<double>(pl, model);
112
113 // Integrate to timeMax
114 bool integratorStatus = integrator->advanceTime();
115 TEST_ASSERT(integratorStatus);
116
117 // Test if at 'Final Time'
118 double time = integrator->getTime();
119 double timeFinal =pl->sublist("Demo Integrator")
120 .sublist("Time Step Control").get<double>("Final Time");
121 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
122
123 // Time-integrated solution and the exact solution using the fact that g = x
124 RCP<const Thyra::VectorBase<double> > x_vec = integrator->getX();
125 RCP<const Thyra::MultiVectorBase<double> > DxDp_vec = integrator->getDgDp();
126 const double x = Thyra::get_ele(*x_vec, 0);
127 const double dxdb = Thyra::get_ele(*(DxDp_vec->col(0)), 0);
128 const double x_exact = model->getSteadyStateSolution();
129 const double dxdb_exact = model->getSteadyStateSolutionSensitivity();
130
131 TEST_FLOATING_EQUALITY( x, x_exact, 1.0e-6 );
132 TEST_FLOATING_EQUALITY( dxdb, dxdb_exact, 1.0e-6 );
133}
134
135
136} // namespace Tempus_Test
Simple quadratic equation with a stable steady-state. This is a simple differential equation.
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
void test_pseudotransient_fsa(const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)