Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_DIRK_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 "Thyra_VectorStdOps.hpp"
16#include "Thyra_MultiVectorStdOps.hpp"
17
18#include "Tempus_IntegratorBasic.hpp"
19#include "Tempus_IntegratorAdjointSensitivity.hpp"
20
21#include "Thyra_DefaultMultiVectorProductVector.hpp"
22#include "Thyra_DefaultProductVector.hpp"
23
24#include "../TestModels/SinCosModel.hpp"
25#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26
27#include <fstream>
28#include <vector>
29
30namespace Tempus_Test {
31
32using Teuchos::getParametersFromXmlFile;
33using Teuchos::ParameterList;
34using Teuchos::RCP;
35using Teuchos::sublist;
36
37using Tempus::IntegratorBasic;
38using Tempus::SolutionHistory;
39using Tempus::SolutionState;
40
41// ************************************************************
42// ************************************************************
43TEUCHOS_UNIT_TEST(DIRK, SinCos_ASA)
44{
45 std::vector<std::string> RKMethods;
46 RKMethods.push_back("General DIRK");
47 RKMethods.push_back("RK Backward Euler");
48 RKMethods.push_back("DIRK 1 Stage Theta Method");
49 RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
50 RKMethods.push_back("SDIRK 2 Stage 2nd order");
51 RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
52 RKMethods.push_back("SDIRK 2 Stage 3rd order");
53 RKMethods.push_back("EDIRK 2 Stage 3rd order");
54 RKMethods.push_back("EDIRK 2 Stage Theta Method");
55 RKMethods.push_back("SDIRK 3 Stage 4th order");
56 RKMethods.push_back("SDIRK 5 Stage 4th order");
57 RKMethods.push_back("SDIRK 5 Stage 5th order");
58
59 std::vector<double> RKMethodErrors;
60 RKMethodErrors.push_back(8.48235e-05);
61 RKMethodErrors.push_back(0.0383339);
62 RKMethodErrors.push_back(0.000221028);
63 RKMethodErrors.push_back(0.0428449);
64 RKMethodErrors.push_back(8.48235e-05);
65 RKMethodErrors.push_back(0.000297933);
66 RKMethodErrors.push_back(4.87848e-06);
67 RKMethodErrors.push_back(7.30827e-07);
68 RKMethodErrors.push_back(0.000272997);
69 RKMethodErrors.push_back(3.10132e-07);
70 RKMethodErrors.push_back(7.56838e-10);
71 RKMethodErrors.push_back(1.32374e-10);
72
73 Teuchos::RCP<const Teuchos::Comm<int> > comm =
74 Teuchos::DefaultComm<int>::getComm();
75
76 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
77 std::string RKMethod_ = RKMethods[m];
78 std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
79 std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
80 std::vector<double> StepSize;
81 std::vector<double> ErrorNorm;
82 const int nTimeStepSizes = 2; // 7 for error plots
83 double dt = 0.05;
84 double order = 0.0;
85 for (int n = 0; n < nTimeStepSizes; n++) {
86 // Read params from .xml file
87 RCP<ParameterList> pList =
88 getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
89
90 // Setup the SinCosModel
91 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
92 RCP<SinCosModel<double> > model =
93 Teuchos::rcp(new SinCosModel<double>(scm_pl));
94
95 // Set the Stepper
96 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
97 pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
98 if (RKMethods[m] == "DIRK 1 Stage Theta Method" ||
99 RKMethods[m] == "EDIRK 2 Stage Theta Method") {
100 pl->sublist("Default Stepper").set<double>("theta", 0.5);
101 }
102 else if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
103 pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
104 }
105 else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
106 pl->sublist("Default Stepper")
107 .set<std::string>("Gamma Type", "3rd Order A-stable");
108 }
109
110 dt /= 2;
111
112 // Setup sensitivities
113 ParameterList& sens_pl = pl->sublist("Sensitivities");
114 sens_pl.set("Mass Matrix Is Identity", true); // Necessary for explicit
115 ParameterList& interp_pl = pl->sublist("Default Integrator")
116 .sublist("Solution History")
117 .sublist("Interpolator");
118 interp_pl.set("Interpolator Type", "Lagrange");
119 interp_pl.set("Order", 4); // All RK methods here are at most 5th order
120
121 // Setup the Integrator and reset initial time step
122 pl->sublist("Default Integrator")
123 .sublist("Time Step Control")
124 .set("Initial Time Step", dt);
125 RCP<Tempus::IntegratorAdjointSensitivity<double> > integrator =
126 Tempus::createIntegratorAdjointSensitivity<double>(pl, model);
127 order = integrator->getStepper()->getOrder();
128
129 // Initial Conditions
130 // During the Integrator construction, the initial SolutionState
131 // is set by default to model->getNominalVales().get_x(). However,
132 // the application can set it also by
133 // integrator->initializeSolutionHistory.
134 RCP<Thyra::VectorBase<double> > x0 =
135 model->getNominalValues().get_x()->clone_v();
136 const int num_param = model->get_p_space(0)->dim();
137 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
138 Thyra::createMembers(model->get_x_space(), num_param);
139 for (int i = 0; i < num_param; ++i)
140 Thyra::assign(DxDp0->col(i).ptr(),
141 *(model->getExactSensSolution(i, 0.0).get_x()));
142 integrator->initializeSolutionHistory(0.0, x0, Teuchos::null,
143 Teuchos::null, DxDp0, Teuchos::null,
144 Teuchos::null);
145
146 // Integrate to timeMax
147 bool integratorStatus = integrator->advanceTime();
148 TEST_ASSERT(integratorStatus)
149
150 // Test if at 'Final Time'
151 double time = integrator->getTime();
152 double timeFinal = pl->sublist("Default Integrator")
153 .sublist("Time Step Control")
154 .get<double>("Final Time");
155 double tol = 100.0 * std::numeric_limits<double>::epsilon();
156 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
157
158 // Time-integrated solution and the exact solution along with
159 // sensitivities (relying on response g(x) = x). Note we must transpose
160 // dg/dp since the integrator returns it in gradient form.
161 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
162 RCP<const Thyra::MultiVectorBase<double> > DgDp = integrator->getDgDp();
163 RCP<Thyra::MultiVectorBase<double> > DxDp =
164 Thyra::createMembers(model->get_x_space(), num_param);
165 {
166 Thyra::ConstDetachedMultiVectorView<double> dgdp_view(*DgDp);
167 Thyra::DetachedMultiVectorView<double> dxdp_view(*DxDp);
168 const int num_g = DgDp->domain()->dim();
169 for (int i = 0; i < num_g; ++i)
170 for (int j = 0; j < num_param; ++j) dxdp_view(i, j) = dgdp_view(j, i);
171 }
172 RCP<const Thyra::VectorBase<double> > x_exact =
173 model->getExactSolution(time).get_x();
174 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
175 Thyra::createMembers(model->get_x_space(), num_param);
176 for (int i = 0; i < num_param; ++i)
177 Thyra::assign(DxDp_exact->col(i).ptr(),
178 *(model->getExactSensSolution(i, time).get_x()));
179
180 // Plot sample solution and exact solution
181 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
182 typedef Thyra::DefaultProductVector<double> DPV;
183 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
184
185 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_AdjSens.dat");
186 RCP<const SolutionHistory<double> > solutionHistory =
187 integrator->getSolutionHistory();
188 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
189 RCP<const SolutionState<double> > solutionState =
190 (*solutionHistory)[i];
191 const double time_i = solutionState->getTime();
192 RCP<const DPV> x_prod_plot =
193 Teuchos::rcp_dynamic_cast<const DPV>(solutionState->getX());
194 RCP<const Thyra::VectorBase<double> > x_plot =
195 x_prod_plot->getVectorBlock(0);
196 RCP<const DMVPV> adjoint_prod_plot =
197 Teuchos::rcp_dynamic_cast<const DMVPV>(
198 x_prod_plot->getVectorBlock(1));
199 RCP<const Thyra::MultiVectorBase<double> > adjoint_plot =
200 adjoint_prod_plot->getMultiVector();
201 RCP<const Thyra::VectorBase<double> > x_exact_plot =
202 model->getExactSolution(time_i).get_x();
203 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
204 << get_ele(*(x_plot), 0) << std::setw(11)
205 << get_ele(*(x_plot), 1) << std::setw(11)
206 << get_ele(*(adjoint_plot->col(0)), 0) << std::setw(11)
207 << get_ele(*(adjoint_plot->col(0)), 1) << std::setw(11)
208 << get_ele(*(adjoint_plot->col(1)), 0) << std::setw(11)
209 << get_ele(*(adjoint_plot->col(1)), 1) << std::setw(11)
210 << get_ele(*(x_exact_plot), 0) << std::setw(11)
211 << get_ele(*(x_exact_plot), 1) << std::endl;
212 }
213 ftmp.close();
214 }
215
216 // Calculate the error
217 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
218 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
219 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
220 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
221 StepSize.push_back(dt);
222 double L2norm = Thyra::norm_2(*xdiff);
223 L2norm *= L2norm;
224 Teuchos::Array<double> L2norm_DxDp(num_param);
225 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
226 for (int i = 0; i < num_param; ++i)
227 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
228 L2norm = std::sqrt(L2norm);
229 ErrorNorm.push_back(L2norm);
230
231 // out << " n = " << n << " dt = " << dt << " error = " << L2norm
232 // << std::endl;
233 }
234
235 if (comm->getRank() == 0) {
236 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_AdjSens-Error.dat");
237 double error0 = 0.8 * ErrorNorm[0];
238 for (int n = 0; n < (int)StepSize.size(); n++) {
239 ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
240 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
241 }
242 ftmp.close();
243 }
244
245 // if (RKMethods[m] == "SDIRK 5 Stage 4th order") {
246 // StepSize.pop_back(); StepSize.pop_back();
247 // ErrorNorm.pop_back(); ErrorNorm.pop_back();
248 // } else if (RKMethods[m] == "SDIRK 5 Stage 5th order") {
249 // StepSize.pop_back(); StepSize.pop_back(); StepSize.pop_back();
250 // ErrorNorm.pop_back(); ErrorNorm.pop_back(); ErrorNorm.pop_back();
251 // }
252
253 // Check the order and intercept
254 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
255 out << " Stepper = " << RKMethods[m] << std::endl;
256 out << " =========================" << std::endl;
257 out << " Expected order: " << order << std::endl;
258 out << " Observed order: " << slope << std::endl;
259 out << " =========================" << std::endl;
260 TEST_FLOATING_EQUALITY(slope, order, 0.03);
261 TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 5.0e-4);
262 }
263 Teuchos::TimeMonitor::summarize();
264}
265
266} // namespace Tempus_Test
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition Tempus_SolutionHistory_decl.hpp:119
Tempus::SolutionState
Solution state for integrators and steppers.
Definition Tempus_SolutionState_decl.hpp:139
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:32
Tempus_Test::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
Definition Tempus_BackwardEuler_ASA.cpp:45
Generated on Thu Oct 9 2025 20:59:09 for Tempus by doxygen 1.9.8