Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_ExplicitRK_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("RK Forward Euler");
50 RKMethods.push_back("RK Explicit 4 Stage");
51 RKMethods.push_back("RK Explicit 3/8 Rule");
52 RKMethods.push_back("RK Explicit 4 Stage 3rd order by Runge");
53 RKMethods.push_back("RK Explicit 5 Stage 3rd order by Kinnmark and Gray");
54 RKMethods.push_back("RK Explicit 3 Stage 3rd order");
55 RKMethods.push_back("RK Explicit 3 Stage 3rd order TVD");
56 RKMethods.push_back("RK Explicit 3 Stage 3rd order by Heun");
57 RKMethods.push_back("RK Explicit Midpoint");
58 RKMethods.push_back("RK Explicit Trapezoidal");
59 RKMethods.push_back("Heuns Method");
60 RKMethods.push_back("General ERK");
61
62 // Check that method_name is valid
63 if (method_name != "") {
64 auto it = std::find(RKMethods.begin(), RKMethods.end(), method_name);
65 TEUCHOS_TEST_FOR_EXCEPTION(
66 it == RKMethods.end(), std::logic_error,
67 "Invalid RK method name '" << method_name << "'");
68 }
69
70 std::vector<double> RKMethodErrors;
71 if (use_combined_method) {
72 RKMethodErrors.push_back(0.183799);
73 RKMethodErrors.push_back(6.88637e-06);
74 RKMethodErrors.push_back(6.88637e-06);
75 RKMethodErrors.push_back(0.000264154);
76 RKMethodErrors.push_back(5.22798e-05);
77 RKMethodErrors.push_back(0.000261896);
78 RKMethodErrors.push_back(0.000261896);
79 RKMethodErrors.push_back(0.000261896);
80 RKMethodErrors.push_back(0.00934377);
81 RKMethodErrors.push_back(0.00934377);
82 RKMethodErrors.push_back(0.00934377);
83 RKMethodErrors.push_back(6.88637e-06);
84 }
85 else {
86 RKMethodErrors.push_back(0.183799);
87 RKMethodErrors.push_back(2.1915e-05);
88 RKMethodErrors.push_back(2.23367e-05);
89 RKMethodErrors.push_back(0.000205051);
90 RKMethodErrors.push_back(2.85141e-05);
91 RKMethodErrors.push_back(0.000126478);
92 RKMethodErrors.push_back(9.64964e-05);
93 RKMethodErrors.push_back(0.000144616);
94 RKMethodErrors.push_back(0.00826159);
95 RKMethodErrors.push_back(0.00710492);
96 RKMethodErrors.push_back(0.00710492);
97 RKMethodErrors.push_back(2.1915e-05);
98 }
99 Teuchos::RCP<const Teuchos::Comm<int> > comm =
100 Teuchos::DefaultComm<int>::getComm();
101
102 for (std::vector<std::string>::size_type m = 0; m != RKMethods.size(); m++) {
103 // If we were given a method to run, skip this method if it doesn't match
104 if (method_name != "" && RKMethods[m] != method_name) continue;
105
106 std::string RKMethod_ = RKMethods[m];
107 std::replace(RKMethod_.begin(), RKMethod_.end(), ' ', '_');
108 std::replace(RKMethod_.begin(), RKMethod_.end(), '/', '.');
109 std::vector<double> StepSize;
110 std::vector<double> ErrorNorm;
111 const int nTimeStepSizes = 7;
112 double dt = 0.2;
113 double order = 0.0;
114 for (int n = 0; n < nTimeStepSizes; n++) {
115 // Read params from .xml file
116 RCP<ParameterList> pList =
117 getParametersFromXmlFile("Tempus_ExplicitRK_SinCos.xml");
118
119 // Setup the SinCosModel
120 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
121 scm_pl->set("Use DfDp as Tangent", use_dfdp_as_tangent);
122 RCP<SinCosModel<double> > model =
123 Teuchos::rcp(new SinCosModel<double>(scm_pl));
124
125 // Set the Stepper
126 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
127 if (RKMethods[m] == "General ERK") {
128 pl->sublist("Demo Integrator").set("Stepper Name", "Demo Stepper 2");
129 }
130 else {
131 pl->sublist("Demo Stepper").set("Stepper Type", RKMethods[m]);
132 }
133
134 dt /= 2;
135
136 // Setup sensitivities
137 ParameterList& sens_pl = pl->sublist("Sensitivities");
138 if (use_combined_method)
139 sens_pl.set("Sensitivity Method", "Combined");
140 else
141 sens_pl.set("Sensitivity Method", "Staggered");
142 sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
143 ParameterList& interp_pl = pl->sublist("Demo Integrator")
144 .sublist("Solution History")
145 .sublist("Interpolator");
146 interp_pl.set("Interpolator Type", "Lagrange");
147 interp_pl.set("Order", 3); // All RK methods here are at most 4th order
148
149 // Setup the Integrator and reset initial time step
150 pl->sublist("Demo Integrator")
151 .sublist("Time Step Control")
152 .set("Initial Time Step", dt);
153 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
154 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
155 order = integrator->getStepper()->getOrder();
156
157 // Initial Conditions
158 double t0 = pl->sublist("Demo Integrator")
159 .sublist("Time Step Control")
160 .get<double>("Initial Time");
161 // RCP<const Thyra::VectorBase<double> > x0 =
162 // model->getExactSolution(t0).get_x()->clone_v();
163 RCP<Thyra::VectorBase<double> > x0 =
164 model->getNominalValues().get_x()->clone_v();
165 const int num_param = model->get_p_space(0)->dim();
166 RCP<Thyra::MultiVectorBase<double> > DxDp0 =
167 Thyra::createMembers(model->get_x_space(), num_param);
168 for (int i = 0; i < num_param; ++i)
169 Thyra::assign(DxDp0->col(i).ptr(),
170 *(model->getExactSensSolution(i, t0).get_x()));
171 integrator->initializeSolutionHistory(t0, x0, Teuchos::null,
172 Teuchos::null, DxDp0, Teuchos::null,
173 Teuchos::null);
174
175 // Integrate to timeMax
176 bool integratorStatus = integrator->advanceTime();
177 TEST_ASSERT(integratorStatus)
178
179 // Test if at 'Final Time'
180 double time = integrator->getTime();
181 double timeFinal = pl->sublist("Demo Integrator")
182 .sublist("Time Step Control")
183 .get<double>("Final Time");
184 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
185
186 // Time-integrated solution and the exact solution
187 RCP<const Thyra::VectorBase<double> > x = integrator->getX();
188 RCP<const Thyra::MultiVectorBase<double> > DxDp = integrator->getDxDp();
189 RCP<const Thyra::VectorBase<double> > x_exact =
190 model->getExactSolution(time).get_x();
191 RCP<Thyra::MultiVectorBase<double> > DxDp_exact =
192 Thyra::createMembers(model->get_x_space(), num_param);
193 for (int i = 0; i < num_param; ++i)
194 Thyra::assign(DxDp_exact->col(i).ptr(),
195 *(model->getExactSensSolution(i, time).get_x()));
196
197 // Plot sample solution and exact solution
198 if (comm->getRank() == 0 && n == nTimeStepSizes - 1) {
199 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
200
201 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens.dat");
202 RCP<const SolutionHistory<double> > solutionHistory =
203 integrator->getSolutionHistory();
204 RCP<Thyra::MultiVectorBase<double> > DxDp_exact_plot =
205 Thyra::createMembers(model->get_x_space(), num_param);
206 for (int i = 0; i < solutionHistory->getNumStates(); i++) {
207 RCP<const SolutionState<double> > solutionState =
208 (*solutionHistory)[i];
209 double time_i = solutionState->getTime();
210 RCP<const DMVPV> x_prod_plot =
211 Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
212 RCP<const Thyra::VectorBase<double> > x_plot =
213 x_prod_plot->getMultiVector()->col(0);
214 RCP<const Thyra::MultiVectorBase<double> > DxDp_plot =
215 x_prod_plot->getMultiVector()->subView(
216 Teuchos::Range1D(1, num_param));
217 RCP<const Thyra::VectorBase<double> > x_exact_plot =
218 model->getExactSolution(time_i).get_x();
219 for (int j = 0; j < num_param; ++j)
220 Thyra::assign(DxDp_exact_plot->col(j).ptr(),
221 *(model->getExactSensSolution(j, time_i).get_x()));
222 ftmp << std::fixed << std::setprecision(7) << time_i << std::setw(11)
223 << get_ele(*(x_plot), 0) << std::setw(11)
224 << get_ele(*(x_plot), 1);
225 for (int j = 0; j < num_param; ++j)
226 ftmp << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 0)
227 << std::setw(11) << get_ele(*(DxDp_plot->col(j)), 1);
228 ftmp << std::setw(11) << get_ele(*(x_exact_plot), 0) << std::setw(11)
229 << get_ele(*(x_exact_plot), 1);
230 for (int j = 0; j < num_param; ++j)
231 ftmp << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 0)
232 << std::setw(11) << get_ele(*(DxDp_exact_plot->col(j)), 1);
233 ftmp << std::endl;
234 }
235 ftmp.close();
236 }
237
238 // Calculate the error
239 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
240 RCP<Thyra::MultiVectorBase<double> > DxDpdiff = DxDp->clone_mv();
241 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
242 Thyra::V_VmV(DxDpdiff.ptr(), *DxDp_exact, *DxDp);
243 StepSize.push_back(dt);
244 double L2norm = Thyra::norm_2(*xdiff);
245 L2norm *= L2norm;
246 Teuchos::Array<double> L2norm_DxDp(num_param);
247 Thyra::norms_2(*DxDpdiff, L2norm_DxDp());
248 for (int i = 0; i < num_param; ++i)
249 L2norm += L2norm_DxDp[i] * L2norm_DxDp[i];
250 L2norm = std::sqrt(L2norm);
251 ErrorNorm.push_back(L2norm);
252
253 out << " n = " << n << " dt = " << dt << " error = " << L2norm
254 << std::endl;
255 }
256
257 // Check the order and intercept
258 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
259 out << " Stepper = " << RKMethods[m] << std::endl;
260 out << " =========================" << std::endl;
261 out << " Expected order: " << order << std::endl;
262 out << " Observed order: " << slope << std::endl;
263 out << " =========================" << std::endl;
264 TEST_FLOATING_EQUALITY(slope, order, 0.04);
265 TEST_FLOATING_EQUALITY(ErrorNorm[0], RKMethodErrors[m], 1.0e-4);
266
267 if (comm->getRank() == 0) {
268 std::ofstream ftmp("Tempus_" + RKMethod_ + "_SinCos_Sens-Error.dat");
269 double error0 = 0.8 * ErrorNorm[0];
270 for (int n = 0; n < nTimeStepSizes; n++) {
271 ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
272 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
273 }
274 ftmp.close();
275 }
276 }
277
278 Teuchos::TimeMonitor::summarize();
279}
280
281} // 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
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