Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_DIRK_FSA.hpp
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_IntegratorForwardSensitivity.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// ************************************************************
43void test_sincos_fsa(const std::string& method_name,
44 const bool use_combined_method,
45 const bool use_dfdp_as_tangent, Teuchos::FancyOStream& out,
46 bool& success)
47{
48 std::vector<std::string> RKMethods;
49 RKMethods.push_back("General DIRK");
50 RKMethods.push_back("RK Backward Euler");
51 RKMethods.push_back("DIRK 1 Stage Theta Method");
52 RKMethods.push_back("RK Implicit 1 Stage 1st order Radau IA");
53 RKMethods.push_back("RK Implicit Midpoint");
54 RKMethods.push_back("SDIRK 2 Stage 2nd order");
55 RKMethods.push_back("RK Implicit 2 Stage 2nd order Lobatto IIIB");
56 RKMethods.push_back("SDIRK 2 Stage 3rd order");
57 RKMethods.push_back("EDIRK 2 Stage 3rd order");
58 RKMethods.push_back("EDIRK 2 Stage Theta Method");
59 RKMethods.push_back("SDIRK 3 Stage 4th order");
60 RKMethods.push_back("SDIRK 5 Stage 4th order");
61 RKMethods.push_back("SDIRK 5 Stage 5th order");
62 RKMethods.push_back("SDIRK 2(1) Pair");
63 RKMethods.push_back("RK Trapezoidal Rule");
64 RKMethods.push_back("RK Crank-Nicolson");
65
66 // Check that method_name is valid
67 if (method_name != "") {
68 auto it = std::find(RKMethods.begin(), RKMethods.end(), method_name);
69 TEUCHOS_TEST_FOR_EXCEPTION(
70 it == RKMethods.end(), std::logic_error,
71 "Invalid RK method name '" << method_name << "'");
72 }
73
74 std::vector<double> RKMethodErrors;
75 if (use_combined_method) {
76 RKMethodErrors.push_back(0.000144507);
77 RKMethodErrors.push_back(0.0428449);
78 RKMethodErrors.push_back(0.000297933);
79 RKMethodErrors.push_back(0.0428449);
80 RKMethodErrors.push_back(0.000297933);
81 RKMethodErrors.push_back(0.000144507);
82 RKMethodErrors.push_back(0.000297933);
83 RKMethodErrors.push_back(8.65434e-06);
84 RKMethodErrors.push_back(1.3468e-06);
85 RKMethodErrors.push_back(0.000297933);
86 RKMethodErrors.push_back(5.44037e-07);
87 RKMethodErrors.push_back(2.77342e-09);
88 RKMethodErrors.push_back(1.21689e-10);
89 RKMethodErrors.push_back(0.000603848);
90 RKMethodErrors.push_back(0.000297933);
91 RKMethodErrors.push_back(0.000297933);
92 }
93 else {
94 RKMethodErrors.push_back(0.000125232);
95 RKMethodErrors.push_back(0.0428449);
96 RKMethodErrors.push_back(0.000221049);
97 RKMethodErrors.push_back(0.0383339);
98 RKMethodErrors.push_back(0.000221049);
99 RKMethodErrors.push_back(0.000125232);
100 RKMethodErrors.push_back(0.000272997);
101 RKMethodErrors.push_back(4.79475e-06);
102 RKMethodErrors.push_back(9.63899e-07);
103 RKMethodErrors.push_back(0.000297933);
104 RKMethodErrors.push_back(2.9362e-07);
105 RKMethodErrors.push_back(9.20081e-08);
106 RKMethodErrors.push_back(9.16252e-08);
107 RKMethodErrors.push_back(0.00043969);
108 RKMethodErrors.push_back(0.000297933);
109 RKMethodErrors.push_back(0.000297933);
110 }
111
112 Teuchos::RCP<const Teuchos::Comm<int> > comm =
113 Teuchos::DefaultComm<int>::getComm();
114
115 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
116 // If we were given a method to run, skip this method if it doesn't match
117 if (method_name != "" && RKMethods[m] != method_name) continue;
118
119 std::string RKMethod_ = RKMethods[m];
120 std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
121 std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
122 std::vector<double> StepSize;
123 std::vector<double> ErrorNorm;
124 const int nTimeStepSizes = 3; // 7 for error plots
125 double dt = 0.05;
126 double order = 0.0;
127 for (int n = 0; n < nTimeStepSizes; n++) {
128 // Read params from .xml file
129 RCP<ParameterList> pList =
130 getParametersFromXmlFile("Tempus_DIRK_SinCos.xml");
131
132 // Setup the SinCosModel
133 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
134 scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
135 RCP<SinCosModel<double> > model =
136 Teuchos::rcp(new SinCosModel<double>(scm_pl));
137
138 // Set the Stepper
139 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
140 pl->sublist("Default Stepper").set("Stepper Type", RKMethods[m]);
141 if (RKMethods[m] == "SDIRK 2 Stage 2nd order") {
142 pl->sublist("Default Stepper").set("gamma", 0.2928932188134524);
143 }
144 else if (RKMethods[m] == "SDIRK 2 Stage 3rd order") {
145 pl->sublist("Default Stepper")
146 .set<std::string>("Gamma Type", "3rd Order A-stable");
147 }
148
149 dt /= 2;
150
151 // Setup sensitivities
152 ParameterList& sens_pl = pl->sublist("Sensitivities");
153 if (use_combined_method)
154 sens_pl.set("Sensitivity Method", "Combined");
155 else
156 sens_pl.set("Sensitivity Method", "Staggered");
157 sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
158 ParameterList& interp_pl = pl->sublist("Default Integrator")
159 .sublist("Solution History")
160 .sublist("Interpolator");
161 interp_pl.set("Interpolator Type", "Lagrange");
162 interp_pl.set("Order", 4); // All RK methods here are at most 5th order
163
164 // Setup the Integrator and reset initial time step
165 pl->sublist("Default Integrator")
166 .sublist("Time Step Control")
167 .set("Initial Time Step", dt);
168 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
169 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
170 order = integrator->getStepper()->getOrder();
171
172 // Initial Conditions
173 // During the Integrator construction, the initial SolutionState
174 // is set by default to model->getNominalVales().get_x(). However,
175 // the application can set it also by
176 // integrator->initializeSolutionHistory.
177 RCP<Thyra::VectorBase<double> > x0 =
178 model->getNominalValues().get_x()->clone_v();
179 const int num_param = model->get_p_space(0)->dim();
180 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
181 Thyra::createMembers(model->get_x_space(), num_param);
182 for (int i = 0; i < num_param; ++i)
183 Thyra::assign(DxDp0->col(i).ptr(),
184 *(model->getExactSensSolution(i, 0.0).get_x()));
185 integrator->initializeSolutionHistory(0.0, x0, Teuchos::null,
186 Teuchos::null, DxDp0, Teuchos::null,
187 Teuchos::null);
188
189 // Integrate to timeMax
190 bool integratorStatus = integrator->advanceTime();
191 TEST_ASSERT(integratorStatus)
192
193 // Test if at 'Final Time'
194 double time = integrator->getTime();
195 double timeFinal = pl->sublist("Default Integrator")
196 .sublist("Time Step Control")
197 .get<double>("Final Time");
198 double tol = 100.0 * std::numeric_limits<double>::epsilon();
199 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
200
201 // Time-integrated solution and the exact solution
202 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
203 RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
204 RCP<const Thyra::VectorBase<double> > x_exact =
205 model->getExactSolution(time).get_x();
206 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
207 Thyra::createMembers(model->get_x_space(), num_param);
208 for (int i = 0; i < num_param; ++i)
209 Thyra::assign(DxDp_exact->col(i).ptr(),
210 *(model->getExactSensSolution(i, time).get_x()));
211
212 // Plot sample solution and exact solution
213 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
214 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
215
216 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens.dat");
217 RCP<const SolutionHistory<double> > solutionHistory =
218 integrator->getSolutionHistory();
219 RCP<Thyra::MultiVectorBase<double> > DxDp_exact_plot =
220 Thyra::createMembers(model->get_x_space(), num_param);
221 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
222 RCP<const SolutionState<double> > solutionState =
223 (*solutionHistory)[i];
224 double time_i = solutionState->getTime();
225 RCP<const DMVPV> x_prod_plot =
226 Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
227 RCP<const Thyra::VectorBase<double> > x_plot =
228 x_prod_plot->getMultiVector()->col(0);
229 RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
230 x_prod_plot->getMultiVector()->subView(
231 Teuchos::Range1D(1, num_param));
232 RCP<const Thyra::VectorBase<double> > x_exact_plot =
233 model->getExactSolution(time_i).get_x();
234 for (int j = 0; j < num_param; ++j)
235 Thyra::assign(DxDp_exact_plot->col(j).ptr(),
236 *(model->getExactSensSolution(j, time_i).get_x()));
237 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
238 << get_ele(*(x_plot), 0) << std::setw(11)
239 << get_ele(*(x_plot), 1);
240 for (int j = 0; j < num_param; ++j)
241 ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
242 << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
243 ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
244 << get_ele(*(x_exact_plot), 1);
245 for (int j = 0; j < num_param; ++j)
246 ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
247 << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
248 ftmp << std::endl;
249 }
250 ftmp.close();
251 }
252
253 // Calculate the error
254 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
255 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
256 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
257 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
258 StepSize.push_back(dt);
259 double L2norm = Thyra::norm_2(*xdiff);
260 L2norm *= L2norm;
261 Teuchos::Array<double> L2norm_DxDp(num_param);
262 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
263 for (int i = 0; i < num_param; ++i)
264 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
265 L2norm = std::sqrt(L2norm);
266 ErrorNorm.push_back(L2norm);
267
268 // out << " n = " << n << " dt = " << dt << " error = " << L2norm
269 // << std::endl;
270 }
271
272 if (comm->getRank() == 0) {
273 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens-Error.dat");
274 double error0 = 0.8 * ErrorNorm[0];
275 for (int n = 0; n < (int)StepSize.size(); n++) {
276 ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
277 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
278 }
279 ftmp.close();
280 }
281
282 // if (RKMethods[m] == "SDIRK 5 Stage 4th order") {
283 // StepSize.pop_back(); StepSize.pop_back();
284 // ErrorNorm.pop_back(); ErrorNorm.pop_back();
285 // } else if (RKMethods[m] == "SDIRK 5 Stage 5th order") {
286 // StepSize.pop_back(); StepSize.pop_back(); StepSize.pop_back();
287 // ErrorNorm.pop_back(); ErrorNorm.pop_back(); ErrorNorm.pop_back();
288 // }
289
290 // Check the order and intercept
291 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
292 out << " Stepper = " << RKMethods[m] << std::endl;
293 out << " =========================" << std::endl;
294 out << " Expected order: " << order << std::endl;
295 out << " Observed order: " << slope << std::endl;
296 out << " =========================" << std::endl;
297
298 // Can only seem to get at most 4th order when using staggered method
299 double order_expected = use_combined_method ? order : std::min(order, 4.0);
300 TEST_FLOATING_EQUALITY(slope, order_expected, 0.03);
301 TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 5.0e-4);
302 }
303 Teuchos::TimeMonitor::summarize();
304}
305
306} // namespace Tempus_Test
method_name
std::string method_name
Definition Tempus_DIRK_Combined_FSA.cpp:16
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::test_sincos_fsa
void test_sincos_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
Definition Tempus_BackwardEuler_FSA.hpp:46
Generated on Thu Oct 9 2025 20:59:09 for Tempus by doxygen 1.9.8