Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IMEX_RK_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#include "Thyra_DefaultMultiVectorProductVector.hpp"
18#include "Thyra_DefaultProductVector.hpp"
19
20#include "Tempus_IntegratorBasic.hpp"
21#include "Tempus_IntegratorForwardSensitivity.hpp"
22#include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
23
24#include "../TestModels/VanDerPolModel.hpp"
25#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
26#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
27#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
28
29#include <fstream>
30#include <vector>
31
32namespace Tempus_Test {
33
34using Teuchos::getParametersFromXmlFile;
35using Teuchos::ParameterList;
36using Teuchos::RCP;
37using Teuchos::sublist;
38
39using Tempus::IntegratorBasic;
40using Tempus::SolutionHistory;
41using Tempus::SolutionState;
42
43// ************************************************************
44// ************************************************************
45void test_vdp_fsa(const bool use_combined_method,
46 const bool use_dfdp_as_tangent, Teuchos::FancyOStream& out,
47 bool& success)
48{
49 std::vector<std::string> stepperTypes;
50 stepperTypes.push_back("IMEX RK 1st order");
51 stepperTypes.push_back("IMEX RK SSP2");
52 stepperTypes.push_back("IMEX RK ARS 233");
53 stepperTypes.push_back("General IMEX RK");
54
55 std::vector<double> stepperOrders;
56 std::vector<double> stepperErrors;
57 if (use_combined_method) {
58 stepperOrders.push_back(1.1198);
59 stepperOrders.push_back(1.98931);
60 stepperOrders.push_back(2.60509);
61 stepperOrders.push_back(1.992);
62
63 stepperErrors.push_back(0.00619674);
64 stepperErrors.push_back(0.294989);
65 stepperErrors.push_back(0.0062125);
66 stepperErrors.push_back(0.142489);
67 }
68 else {
69 stepperOrders.push_back(1.1198);
70 stepperOrders.push_back(2.05232);
71 stepperOrders.push_back(2.43013);
72 stepperOrders.push_back(2.07975);
73
74 stepperErrors.push_back(0.00619674);
75 stepperErrors.push_back(0.0534172);
76 stepperErrors.push_back(0.00224533);
77 stepperErrors.push_back(0.032632);
78 }
79 std::vector<double> stepperInitDt;
80 stepperInitDt.push_back(0.0125);
81 stepperInitDt.push_back(0.05);
82 stepperInitDt.push_back(0.05);
83 stepperInitDt.push_back(0.05);
84
85 Teuchos::RCP<const Teuchos::Comm<int>> comm =
86 Teuchos::DefaultComm<int>::getComm();
87
88 std::vector<std::string>::size_type m;
89 for (m = 0; m != stepperTypes.size(); m++) {
90 std::string stepperType = stepperTypes[m];
91 std::string stepperName = stepperTypes[m];
92 std::replace(stepperName.begin(), stepperName.end(), ' ', '_');
93 std::replace(stepperName.begin(), stepperName.end(), '/', '.');
94
95 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
96 std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
97 std::vector<double> StepSize;
98 std::vector<double> ErrorNorm;
99 const int nTimeStepSizes = 3; // 6 for error plot
100 double dt = stepperInitDt[m];
101 double order = 0.0;
102 for (int n = 0; n < nTimeStepSizes; n++) {
103 // Read params from .xml file
104 RCP<ParameterList> pList =
105 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
106
107 // Setup the explicit VanDerPol ModelEvaluator
108 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
109 vdpmPL->set("Use DfDp as Tangent", use_dfdp_as_tangent);
110 RCP<VanDerPol_IMEX_ExplicitModel<double>> explicitModel =
111 Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
112
113 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
114 RCP<VanDerPol_IMEX_ImplicitModel<double>> implicitModel =
115 Teuchos::rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
116
117 // Setup the IMEX Pair ModelEvaluator
118 RCP<Tempus::WrapperModelEvaluatorPairIMEX_Basic<double>> model =
119 Teuchos::rcp(new Tempus::WrapperModelEvaluatorPairIMEX_Basic<double>(
120 explicitModel, implicitModel));
121
122 // Setup sensitivities
123 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
124 ParameterList& sens_pl = pl->sublist("Sensitivities");
125 if (use_combined_method)
126 sens_pl.set("Sensitivity Method", "Combined");
127 else {
128 sens_pl.set("Sensitivity Method", "Staggered");
129 sens_pl.set("Reuse State Linear Solver", true);
130 }
131 sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
132 ParameterList& interp_pl = pl->sublist("Default Integrator")
133 .sublist("Solution History")
134 .sublist("Interpolator");
135 interp_pl.set("Interpolator Type", "Lagrange");
136 interp_pl.set("Order", 2); // All RK methods here are at most 3rd order
137
138 // Set the Stepper
139 if (stepperType == "General IMEX RK") {
140 // use the appropriate stepper sublist
141 pl->sublist("Default Integrator")
142 .set("Stepper Name", "General IMEX RK");
143 }
144 else {
145 pl->sublist("Default Stepper").set("Stepper Type", stepperType);
146 }
147
148 // Set the step size
149 if (n == nTimeStepSizes - 1)
150 dt /= 10.0;
151 else
152 dt /= 2;
153
154 // Setup the Integrator and reset initial time step
155 pl->sublist("Default Integrator")
156 .sublist("Time Step Control")
157 .set("Initial Time Step", dt);
158 pl->sublist("Default Integrator")
159 .sublist("Time Step Control")
160 .remove("Time Step Control Strategy");
161 RCP<Tempus::IntegratorForwardSensitivity<double>> integrator =
162 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
163 order = integrator->getStepper()->getOrder();
164
165 // Integrate to timeMax
166 bool integratorStatus = integrator->advanceTime();
167 TEST_ASSERT(integratorStatus)
168
169 // Test if at 'Final Time'
170 double time = integrator->getTime();
171 double timeFinal = pl->sublist("Default Integrator")
172 .sublist("Time Step Control")
173 .get<double>("Final Time");
174 double tol = 100.0 * std::numeric_limits<double>::epsilon();
175 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
176
177 // Store off the final solution and step size
178 auto solution = Thyra::createMember(model->get_x_space());
179 auto sensitivity = Thyra::createMember(model->get_x_space());
180 Thyra::copy(*(integrator->getX()), solution.ptr());
181 Thyra::copy(*(integrator->getDxDp()->col(0)), sensitivity.ptr());
182 solutions.push_back(solution);
183 sensitivities.push_back(sensitivity);
184 StepSize.push_back(dt);
185
186 // Output finest temporal solution for plotting
187 if (comm->getRank() == 0 && ((n == 0) || (n == nTimeStepSizes - 1))) {
188 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
189
190 std::string fname = "Tempus_" + stepperName + "_VanDerPol_Sens-Ref.dat";
191 if (n == 0) fname = "Tempus_" + stepperName + "_VanDerPol_Sens.dat";
192 std::ofstream ftmp(fname);
193 RCP<const SolutionHistory<double>> solutionHistory =
194 integrator->getSolutionHistory();
195 int nStates = solutionHistory->getNumStates();
196 for (int i = 0; i < nStates; i++) {
197 RCP<const SolutionState<double>> solutionState =
198 (*solutionHistory)[i];
199 RCP<const DMVPV> x_prod =
200 Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
201 RCP<const Thyra::VectorBase<double>> x =
202 x_prod->getMultiVector()->col(0);
203 RCP<const Thyra::VectorBase<double>> dxdp =
204 x_prod->getMultiVector()->col(1);
205 double ttime = solutionState->getTime();
206 ftmp << std::fixed << std::setprecision(7) << ttime << " "
207 << std::setw(11) << get_ele(*x, 0) << " " << std::setw(11)
208 << get_ele(*x, 1) << " " << std::setw(11) << get_ele(*dxdp, 0)
209 << " " << std::setw(11) << get_ele(*dxdp, 1) << std::endl;
210 }
211 ftmp.close();
212 }
213 }
214
215 // Calculate the error - use the most temporally refined mesh for
216 // the reference solution.
217 auto ref_solution = solutions[solutions.size() - 1];
218 auto ref_sensitivity = sensitivities[solutions.size() - 1];
219 std::vector<double> StepSizeCheck;
220 for (std::size_t i = 0; i < (solutions.size() - 1); ++i) {
221 auto sol = solutions[i];
222 auto sen = sensitivities[i];
223 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
224 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
225 const double L2norm_sol = Thyra::norm_2(*sol);
226 const double L2norm_sen = Thyra::norm_2(*sen);
227 const double L2norm =
228 std::sqrt(L2norm_sol * L2norm_sol + L2norm_sen * L2norm_sen);
229 StepSizeCheck.push_back(StepSize[i]);
230 ErrorNorm.push_back(L2norm);
231
232 out << " n = " << i << " dt = " << StepSize[i] << " error = " << L2norm
233 << std::endl;
234 }
235
236 // Check the order and intercept
237 double slope =
238 computeLinearRegressionLogLog<double>(StepSizeCheck, ErrorNorm);
239 out << " Stepper = " << stepperType << std::endl;
240 out << " =========================" << std::endl;
241 out << " Expected order: " << order << std::endl;
242 out << " Observed order: " << slope << std::endl;
243 out << " =========================" << std::endl;
244 TEST_FLOATING_EQUALITY(slope, stepperOrders[m], 0.02);
245 TEST_FLOATING_EQUALITY(ErrorNorm[0], stepperErrors[m], 1.0e-4);
246
247 // Write error data
248 {
249 std::ofstream ftmp("Tempus_" + stepperName + "_VanDerPol_Sens-Error.dat");
250 double error0 = 0.8 * ErrorNorm[0];
251 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
252 ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
253 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
254 }
255 ftmp.close();
256 }
257 }
258 Teuchos::TimeMonitor::summarize();
259}
260
261} // 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::WrapperModelEvaluatorPairIMEX_Basic
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Definition Tempus_WrapperModelEvaluatorPairIMEX_Basic_decl.hpp:40
Tempus_Test::VanDerPol_IMEX_ExplicitModel
van der Pol model formulated for IMEX.
Definition VanDerPol_IMEX_ExplicitModel_decl.hpp:110
Tempus_Test::VanDerPol_IMEX_ImplicitModel
van der Pol model formulated for IMEX-RK.
Definition VanDerPol_IMEX_ImplicitModel_decl.hpp:110
Tempus_Test
Definition Tempus_BackwardEuler_ASA.cpp:32
Tempus_Test::test_vdp_fsa
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
Definition Tempus_IMEX_RK_FSA.hpp:45
Generated on Thu Oct 9 2025 20:59:09 for Tempus by doxygen 1.9.8