Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_BDF2_ASA.cpp
Go to the documentation of this file.
1//@HEADER
2// *****************************************************************************
3// Tempus: Time Integration and Sensitivity Analysis Package
4//
5// Copyright 2017 NTESS and the Tempus contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8//@HEADER
9
10#include "Teuchos_UnitTestHarness.hpp"
11#include "Teuchos_XMLParameterListHelpers.hpp"
12#include "Teuchos_TimeMonitor.hpp"
13#include "Teuchos_DefaultComm.hpp"
14
15#include "Tempus_config.hpp"
16#include "Tempus_IntegratorBasic.hpp"
17#include "Tempus_IntegratorAdjointSensitivity.hpp"
18
19#include "Thyra_VectorStdOps.hpp"
20#include "Thyra_MultiVectorStdOps.hpp"
21
22#include "../TestModels/SinCosModel.hpp"
23#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
24#include "../TestUtils/Tempus_AdjointResponseTestUtils.hpp"
25
26#include "Thyra_DefaultMultiVectorProductVector.hpp"
27
28#include <fstream>
29#include <limits>
30#include <sstream>
31#include <vector>
32
33namespace Tempus_Test {
34
35using Teuchos::getParametersFromXmlFile;
36using Teuchos::ParameterList;
37using Teuchos::RCP;
38using Teuchos::sublist;
39
40using Tempus::IntegratorBasic;
41using Tempus::SolutionHistory;
42using Tempus::SolutionState;
43
44// ************************************************************
45// ************************************************************
46TEUCHOS_UNIT_TEST(BDF2, SinCos_ASA)
47{
48 std::vector<double> StepSize;
49 std::vector<double> ErrorNorm;
50 const int nTimeStepSizes = 7;
51 double dt = 0.2;
52 double order = 0.0;
53 Teuchos::RCP<const Teuchos::Comm<int> > comm =
54 Teuchos::DefaultComm<int>::getComm();
55 for (int n = 0; n < nTimeStepSizes; n++) {
56 // Read params from .xml file
57 RCP<ParameterList> pList =
58 getParametersFromXmlFile("Tempus_BDF2_SinCos_SA.xml");
59
60 // Setup the SinCosModel
61 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
62 RCP<SinCosModel<double> > model =
63 Teuchos::rcp(new SinCosModel<double>(scm_pl));
64
65 dt /= 2;
66
67 // Setup sensitivities
68 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
69 // ParameterList& sens_pl = pl->sublist("Sensitivities");
70 ParameterList& interp_pl = pl->sublist("Default Integrator")
71 .sublist("Solution History")
72 .sublist("Interpolator");
73 interp_pl.set("Interpolator Type", "Lagrange");
74 interp_pl.set("Order", 1);
75
76 // Setup the Integrator and reset initial time step
77 pl->sublist("Default Integrator")
78 .sublist("Time Step Control")
79 .set("Initial Time Step", dt);
80 RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
81 Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
82 order = integrator->getStepper()->getOrder();
83
84 // Initial Conditions
85 double t0 = pl->sublist("Default Integrator")
86 .sublist("Time Step Control")
87 .get<double>("Initial Time");
88 RCP<const Thyra::VectorBase<double> > x0 =
89 model->getExactSolution(t0).get_x();
90 const int num_param = model->get_p_space(0)->dim();
91 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
92 Thyra::createMembers(model->get_x_space(), num_param);
93 for (int i = 0; i < num_param; ++i)
94 Thyra::assign(DxDp0->col(i).ptr(),
95 *(model->getExactSensSolution(i, t0).get_x()));
96 integrator->initializeSolutionHistory(t0, x0, Teuchos::null, Teuchos::null,
97 DxDp0, Teuchos::null, Teuchos::null);
98
99 // Integrate to timeMax
100 bool integratorStatus = integrator->advanceTime();
101 TEST_ASSERT(integratorStatus)
102
103 // Test if at 'Final Time'
104 double time = integrator->getTime();
105 double timeFinal = pl->sublist("Default Integrator")
106 .sublist("Time Step Control")
107 .get<double>("Final Time");
108 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
109
110 // Time-integrated solution and the exact solution along with
111 // sensitivities (relying on response g(x) = x). Note we must transpose
112 // dg/dp since the integrator returns it in gradient form.
113 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
114 RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
115 RCP<Thyra::MultiVectorBase<double> > DxDp =
116 Thyra::createMembers(model->get_x_space(), num_param);
117 {
118 Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
119 Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
120 const int num_g = DgDp->domain()->dim();
121 for (int i = 0; i < num_g; ++i)
122 for (int j = 0; j < num_param; ++j) dxdp_view(i, j) = dgdp_view(j, i);
123 }
124 RCP<const Thyra::VectorBase<double> > x_exact =
125 model->getExactSolution(time).get_x();
126 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
127 Thyra::createMembers(model->get_x_space(), num_param);
128 for (int i = 0; i < num_param; ++i)
129 Thyra::assign(DxDp_exact->col(i).ptr(),
130 *(model->getExactSensSolution(i, time).get_x()));
131
132 // Plot sample solution and exact solution
133 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
134 typedef Thyra::DefaultProductVector<double> DPV;
135 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
136
137 std::ofstream ftmp("Tempus_BDF2_SinCos_AdjSens.dat");
138 RCP<const SolutionHistory<double> > solutionHistory =
139 integrator->getSolutionHistory();
140 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
141 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
142 const double time_i = solutionState->getTime();
143 RCP<const DPV> x_prod_plot =
144 Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
145 RCP<const Thyra::VectorBase<double> > x_plot =
146 x_prod_plot->getVectorBlock(0);
147 RCP<const DMVPV> adjoint_prod_plot =
148 Teuchos::rcp_dynamic_cast<const DMVPV>(
149 x_prod_plot->getVectorBlock(1));
150 RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
151 adjoint_prod_plot->getMultiVector();
152 RCP<const Thyra::VectorBase<double> > x_exact_plot =
153 model->getExactSolution(time_i).get_x();
154 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
155 << get_ele(*(x_plot), 0) << std::setw(11) << get_ele(*(x_plot), 1)
156 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 0)
157 << std::setw(11) << get_ele(*(adjoint_plot->col(0)), 1)
158 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 0)
159 << std::setw(11) << get_ele(*(adjoint_plot->col(1)), 1)
160 << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
161 << get_ele(*(x_exact_plot), 1) << std::endl;
162 }
163 ftmp.close();
164 }
165
166 // Calculate the error
167 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
168 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
169 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
170 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
171 StepSize.push_back(dt);
172 double L2norm = Thyra::norm_2(*xdiff);
173 L2norm *= L2norm;
174 Teuchos::Array<double> L2norm_DxDp(num_param);
175 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
176 for (int i = 0; i < num_param; ++i)
177 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
178 L2norm = std::sqrt(L2norm);
179 ErrorNorm.push_back(L2norm);
180
181 // out << " n = " << n << " dt = " << dt << " error = " << L2norm
182 // << std::endl;
183 }
184
185 // Check the order and intercept
186 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
187 out << " Stepper = BDF2" << std::endl;
188 out << " =========================" << std::endl;
189 out << " Expected order: " << order << std::endl;
190 out << " Observed order: " << slope << std::endl;
191 out << " =========================" << std::endl;
192 TEST_FLOATING_EQUALITY(slope, 1.95006, 0.015);
193 TEST_FLOATING_EQUALITY(ErrorNorm[0], 0.0137394, 1.0e-4);
194
195 if (comm->getRank() == 0) {
196 std::ofstream ftmp("Tempus_BDF2_SinCos_AdjSens-Error.dat");
197 double error0 = 0.8 * ErrorNorm[0];
198 for (int n = 0; n < nTimeStepSizes; n++) {
199 ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
200 << error0 * (StepSize[n] / StepSize[0]) << std::endl;
201 }
202 ftmp.close();
203 }
204}
205
206// ************************************************************
207// ************************************************************
208// Adjoint sensitivities of the time-integrated response
209// J = int_{t0}^{T} x(t) dt, with BDF2.
210TEUCHOS_UNIT_TEST(BDF2, SinCos_ASA_TimeIntegral)
211{
212 std::vector<double> StepSize, ErrorNorm;
213 double order = 0.0;
214 computeSinCosAdjointResponseErrors(
215 "Tempus_BDF2_SinCos_SA.xml",
216 [](ParameterList& pl) {
217 ParameterList& interp_pl = pl.sublist("Default Integrator")
218 .sublist("Solution History")
219 .sublist("Interpolator");
220 interp_pl.set("Interpolator Type", "Lagrange");
221 interp_pl.set("Order", 1);
222 },
223 "Time Integral", std::vector<double>(), 5, 0.025, false, StepSize, ErrorNorm,
224 order, out, success);
225
226 const double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
227 out << " Expected order: " << order << " Observed order: " << slope << std::endl;
228 TEST_FLOATING_EQUALITY(slope, order, 0.04);
229}
230
231// ************************************************************
232// ************************************************************
233// Adjoint sensitivities of the response J = sum_k x(t_k), including the
234// initial and final times, with BDF2. The adjoint integration is restarted
235// at each time t_k (with the start up stepper of BDF2), so smaller time steps
236// are needed to reach the asymptotic order.
237TEUCHOS_UNIT_TEST(BDF2, SinCos_ASA_SumOverTimes)
238{
239 std::vector<double> StepSize, ErrorNorm;
240 double order = 0.0;
241 computeSinCosAdjointResponseErrors(
242 "Tempus_BDF2_SinCos_SA.xml",
243 [](ParameterList& pl) {
244 ParameterList& interp_pl = pl.sublist("Default Integrator")
245 .sublist("Solution History")
246 .sublist("Interpolator");
247 interp_pl.set("Interpolator Type", "Lagrange");
248 interp_pl.set("Order", 1);
249 },
250 "Sum Over Times", {0.0, 0.3, 0.6, 1.0}, 5, 0.025, false, StepSize, ErrorNorm,
251 order, out, success);
252
253 const double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
254 out << " Expected order: " << order << " Observed order: " << slope << std::endl;
255 TEST_FLOATING_EQUALITY(slope, order, 0.04);
256}
257
258} // namespace Tempus_Test
Tempus_AdjointResponseTestUtils.hpp
Tempus_ConvergenceTestUtils.hpp
Tempus_Test::SinCosModel
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation.
Definition SinCosModel_decl.hpp:110
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:47
Tempus_Test::computeSinCosAdjointResponseErrors
void computeSinCosAdjointResponseErrors(const std::string &xmlFile, const std::function< void(Teuchos::ParameterList &)> &configure, const std::string &responseType, const std::vector< double > &responseTimes, const int nTimeStepSizes, const double dt0, const bool useAdjointModel, std::vector< double > &stepSizes, std::vector< double > &errors, double &stepperOrder, Teuchos::FancyOStream &out, bool &success)
Solve the SinCos adjoint problem with the "Time Integral" or "Sum Over Times" response for a sequence...
Definition Tempus_AdjointResponseTestUtils.hpp:93
Tempus_Test::computeLinearRegressionLogLog
Scalar computeLinearRegressionLogLog(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition Tempus_ConvergenceTestUtils.hpp:159
Tempus::createIntegratorAdjointSensitivity
Teuchos::RCP< IntegratorAdjointSensitivity< Scalar > > createIntegratorAdjointSensitivity()
Nonmember constructor.
Definition Tempus_IntegratorAdjointSensitivity_impl.hpp:900
Generated by doxygen 1.15.0