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: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12#include "Teuchos_DefaultComm.hpp"
13
14#include "Thyra_VectorStdOps.hpp"
15#include "Thyra_MultiVectorStdOps.hpp"
16#include "Thyra_DefaultMultiVectorProductVector.hpp"
17#include "Thyra_DefaultProductVector.hpp"
18
19#include "Tempus_IntegratorBasic.hpp"
20#include "Tempus_IntegratorForwardSensitivity.hpp"
21#include "Tempus_WrapperModelEvaluatorPairIMEX_Basic.hpp"
22
23#include "../TestModels/VanDerPolModel.hpp"
24#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
25#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
27
28#include <fstream>
29#include <vector>
30
31namespace Tempus_Test {
32
33using Teuchos::RCP;
34using Teuchos::ParameterList;
35using Teuchos::sublist;
36using Teuchos::getParametersFromXmlFile;
37
41
42
43// ************************************************************
44// ************************************************************
45void test_vdp_fsa(const bool use_combined_method,
46 const bool use_dfdp_as_tangent,
47 Teuchos::FancyOStream &out, 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 Teuchos::RCP<Teuchos::FancyOStream> my_out =
88 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
89 my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
90 my_out->setOutputToRootOnly(0);
91
92 std::vector<std::string>::size_type m;
93 for(m = 0; m != stepperTypes.size(); m++) {
94
95 std::string stepperType = stepperTypes[m];
96 std::string stepperName = stepperTypes[m];
97 std::replace(stepperName.begin(), stepperName.end(), ' ', '_');
98 std::replace(stepperName.begin(), stepperName.end(), '/', '.');
99
100 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
101 std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
102 std::vector<double> StepSize;
103 std::vector<double> ErrorNorm;
104 const int nTimeStepSizes = 3; // 6 for error plot
105 double dt = stepperInitDt[m];
106 double order = 0.0;
107 for (int n=0; n<nTimeStepSizes; n++) {
108
109 // Read params from .xml file
110 RCP<ParameterList> pList =
111 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
112
113 // Setup the explicit VanDerPol ModelEvaluator
114 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
115 vdpmPL->set("Use DfDp as Tangent", use_dfdp_as_tangent);
116 RCP<VanDerPol_IMEX_ExplicitModel<double> > explicitModel =
117 Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL));
118
119 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
120 RCP<VanDerPol_IMEX_ImplicitModel<double> > implicitModel =
121 Teuchos::rcp(new VanDerPol_IMEX_ImplicitModel<double>(vdpmPL));
122
123 // Setup the IMEX Pair ModelEvaluator
124 RCP<Tempus::WrapperModelEvaluatorPairIMEX_Basic<double> > model =
126 explicitModel, implicitModel));
127
128 // Setup sensitivities
129 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
130 ParameterList& sens_pl = pl->sublist("Sensitivities");
131 if (use_combined_method)
132 sens_pl.set("Sensitivity Method", "Combined");
133 else {
134 sens_pl.set("Sensitivity Method", "Staggered");
135 sens_pl.set("Reuse State Linear Solver", true);
136 }
137 sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
138 ParameterList& interp_pl =
139 pl->sublist("Default Integrator").sublist("Solution History").sublist("Interpolator");
140 interp_pl.set("Interpolator Type", "Lagrange");
141 interp_pl.set("Order", 2); // All RK methods here are at most 3rd order
142
143 // Set the Stepper
144 if (stepperType == "General IMEX RK"){
145 // use the appropriate stepper sublist
146 pl->sublist("Default Integrator").set("Stepper Name", "General IMEX RK");
147 } else {
148 pl->sublist("Default Stepper").set("Stepper Type", stepperType);
149 }
150
151 // Set the step size
152 if (n == nTimeStepSizes-1) dt /= 10.0;
153 else dt /= 2;
154
155 // Setup the Integrator and reset initial time step
156 pl->sublist("Default Integrator")
157 .sublist("Time Step Control").set("Initial Time Step", dt);
158 pl->sublist("Default Integrator")
159 .sublist("Time Step Control").remove("Time Step Control Strategy");
160 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
161 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
162 order = integrator->getStepper()->getOrder();
163
164 // Integrate to timeMax
165 bool integratorStatus = integrator->advanceTime();
166 TEST_ASSERT(integratorStatus)
167
168 // Test if at 'Final Time'
169 double time = integrator->getTime();
170 double timeFinal =pl->sublist("Default Integrator")
171 .sublist("Time Step Control").get<double>("Final Time");
172 double tol = 100.0 * std::numeric_limits<double>::epsilon();
173 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
174
175 // Store off the final solution and step size
176 auto solution = Thyra::createMember(model->get_x_space());
177 auto sensitivity = Thyra::createMember(model->get_x_space());
178 Thyra::copy(*(integrator->getX()),solution.ptr());
179 Thyra::copy(*(integrator->getDxDp()->col(0)),sensitivity.ptr());
180 solutions.push_back(solution);
181 sensitivities.push_back(sensitivity);
182 StepSize.push_back(dt);
183
184 // Output finest temporal solution for plotting
185 if (comm->getRank() == 0 && ((n == 0) || (n == nTimeStepSizes-1))) {
186 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
187
188 std::string fname = "Tempus_"+stepperName+"_VanDerPol_Sens-Ref.dat";
189 if (n == 0) fname = "Tempus_"+stepperName+"_VanDerPol_Sens.dat";
190 std::ofstream ftmp(fname);
191 RCP<const SolutionHistory<double> > solutionHistory =
192 integrator->getSolutionHistory();
193 int nStates = solutionHistory->getNumStates();
194 for (int i=0; i<nStates; i++) {
195 RCP<const SolutionState<double> > solutionState =
196 (*solutionHistory)[i];
197 RCP<const DMVPV> x_prod =
198 Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
199 RCP<const Thyra::VectorBase<double> > x =
200 x_prod->getMultiVector()->col(0);
201 RCP<const Thyra::VectorBase<double> > dxdp =
202 x_prod->getMultiVector()->col(1);
203 double ttime = solutionState->getTime();
204 ftmp << std::fixed << std::setprecision(7)
205 << ttime << " "
206 << std::setw(11) << get_ele(*x, 0) << " "
207 << std::setw(11) << get_ele(*x, 1) << " "
208 << std::setw(11) << get_ele(*dxdp, 0) << " "
209 << std::setw(11) << get_ele(*dxdp, 1)
210 << std::endl;
211 }
212 ftmp.close();
213 }
214 }
215
216 // Calculate the error - use the most temporally refined mesh for
217 // the reference solution.
218 auto ref_solution = solutions[solutions.size()-1];
219 auto ref_sensitivity = sensitivities[solutions.size()-1];
220 std::vector<double> StepSizeCheck;
221 for (std::size_t i=0; i < (solutions.size()-1); ++i) {
222 auto sol = solutions[i];
223 auto sen = sensitivities[i];
224 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
225 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
226 const double L2norm_sol = Thyra::norm_2(*sol);
227 const double L2norm_sen = Thyra::norm_2(*sen);
228 const double L2norm =
229 std::sqrt(L2norm_sol*L2norm_sol + L2norm_sen*L2norm_sen);
230 StepSizeCheck.push_back(StepSize[i]);
231 ErrorNorm.push_back(L2norm);
232
233 *my_out << " n = " << i << " dt = " << StepSize[i]
234 << " error = " << L2norm << std::endl;
235 }
236
237 // Check the order and intercept
238 double slope = 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
262} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)