Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_IMEX_RK_Partitioned_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_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
22
23#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
24#include "../TestModels/VanDerPol_IMEXPart_ImplicitModel.hpp"
26
27#include <fstream>
28#include <vector>
29
30namespace Tempus_Test {
31
32using Teuchos::RCP;
33using Teuchos::ParameterList;
34using Teuchos::sublist;
35using Teuchos::getParametersFromXmlFile;
36
40
41
42// ************************************************************
43// ************************************************************
44void test_vdp_fsa(const std::string& method_name,
45 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("Partitioned IMEX RK 1st order");
51 stepperTypes.push_back("Partitioned IMEX RK SSP2" );
52 stepperTypes.push_back("Partitioned IMEX RK ARS 233" );
53 stepperTypes.push_back("General Partitioned IMEX RK" );
54
55 // Check that method_name is valid
56 if (method_name != "") {
57 auto it = std::find(stepperTypes.begin(), stepperTypes.end(), method_name);
58 TEUCHOS_TEST_FOR_EXCEPTION(it == stepperTypes.end(), std::logic_error,
59 "Invalid stepper type " << method_name);
60 }
61
62 std::vector<double> stepperOrders;
63 std::vector<double> stepperErrors;
64 if (use_dfdp_as_tangent) {
65 if (use_combined_method) {
66 stepperOrders.push_back(1.16082);
67 stepperOrders.push_back(1.97231);
68 stepperOrders.push_back(2.5914);
69 stepperOrders.push_back(1.99148);
70
71 stepperErrors.push_back(0.00820931);
72 stepperErrors.push_back(0.287112);
73 stepperErrors.push_back(0.00646096);
74 stepperErrors.push_back(0.148848);
75 }
76 else {
77 stepperOrders.push_back(1.07932);
78 stepperOrders.push_back(1.97396);
79 stepperOrders.push_back(2.63724);
80 stepperOrders.push_back(1.99133);
81
82 stepperErrors.push_back(0.055626);
83 stepperErrors.push_back(0.198898);
84 stepperErrors.push_back(0.00614135);
85 stepperErrors.push_back(0.0999881);
86 }
87 }
88 else {
89 if (use_combined_method) {
90 stepperOrders.push_back(1.1198);
91 stepperOrders.push_back(1.98931);
92 stepperOrders.push_back(2.60509);
93 stepperOrders.push_back(1.992);
94
95 stepperErrors.push_back(0.00619674);
96 stepperErrors.push_back(0.294989);
97 stepperErrors.push_back(0.0062125);
98 stepperErrors.push_back(0.142489);
99 }
100 else {
101 stepperOrders.push_back(1.07932);
102 stepperOrders.push_back(1.97396);
103 stepperOrders.push_back(2.63724);
104 stepperOrders.push_back(1.99133);
105
106 stepperErrors.push_back(0.055626);
107 stepperErrors.push_back(0.198898);
108 stepperErrors.push_back(0.00614135);
109 stepperErrors.push_back(0.0999881);
110 }
111 }
112
113 std::vector<double> stepperInitDt;
114 stepperInitDt.push_back(0.0125);
115 stepperInitDt.push_back(0.05);
116 stepperInitDt.push_back(0.05);
117 stepperInitDt.push_back(0.05);
118
119 Teuchos::RCP<const Teuchos::Comm<int> > comm =
120 Teuchos::DefaultComm<int>::getComm();
121 Teuchos::RCP<Teuchos::FancyOStream> my_out =
122 Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
123 my_out->setProcRankAndSize(comm->getRank(), comm->getSize());
124 my_out->setOutputToRootOnly(0);
125
126 std::vector<std::string>::size_type m;
127 for(m = 0; m != stepperTypes.size(); m++) {
128
129 // If we were given a method to run, skip this method if it doesn't match
130 if (method_name != "" && stepperTypes[m] != method_name)
131 continue;
132
133 std::string stepperType = stepperTypes[m];
134 std::string stepperName = stepperTypes[m];
135 std::replace(stepperName.begin(), stepperName.end(), ' ', '_');
136 std::replace(stepperName.begin(), stepperName.end(), '/', '.');
137
138 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
139 std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
140 std::vector<double> StepSize;
141 std::vector<double> ErrorNorm;
142 const int nTimeStepSizes = 3; // 6 for error plot
143 double dt = stepperInitDt[m];
144 double order = 0.0;
145 for (int n=0; n<nTimeStepSizes; n++) {
146
147 // Read params from .xml file
148 RCP<ParameterList> pList =
149 getParametersFromXmlFile("Tempus_IMEX_RK_VanDerPol.xml");
150
151 // Setup the explicit VanDerPol ModelEvaluator
152 RCP<ParameterList> vdpmPL = sublist(pList, "VanDerPolModel", true);
153 vdpmPL->set("Use DfDp as Tangent", use_dfdp_as_tangent);
154 const bool useProductVector = true;
155 RCP<VanDerPol_IMEX_ExplicitModel<double> > explicitModel =
156 Teuchos::rcp(new VanDerPol_IMEX_ExplicitModel<double>(vdpmPL,
157 useProductVector));
158
159 // Setup the implicit VanDerPol ModelEvaluator (reuse vdpmPL)
160 RCP<VanDerPol_IMEXPart_ImplicitModel<double> > implicitModel =
161 Teuchos::rcp(new VanDerPol_IMEXPart_ImplicitModel<double>(vdpmPL));
162
163 // Setup the IMEX Pair ModelEvaluator
164 const int numExplicitBlocks = 1;
165 const int parameterIndex = 4;
166 RCP<Tempus::WrapperModelEvaluatorPairPartIMEX_Basic<double> > model =
167 Teuchos::rcp(
169 explicitModel, implicitModel,
170 numExplicitBlocks, parameterIndex));
171
172 // Setup sensitivities
173 RCP<ParameterList> pl = sublist(pList, "Tempus", true);
174 ParameterList& sens_pl = pl->sublist("Sensitivities");
175 if (use_combined_method)
176 sens_pl.set("Sensitivity Method", "Combined");
177 else {
178 sens_pl.set("Sensitivity Method", "Staggered");
179 sens_pl.set("Reuse State Linear Solver", true);
180 }
181 sens_pl.set("Use DfDp as Tangent", use_dfdp_as_tangent);
182 ParameterList& interp_pl =
183 pl->sublist("Default Integrator").sublist("Solution History").sublist("Interpolator");
184 interp_pl.set("Interpolator Type", "Lagrange");
185 interp_pl.set("Order", 2); // All RK methods here are at most 3rd order
186
187 // Set the Stepper
188 if (stepperType == "General Partitioned IMEX RK"){
189 // use the appropriate stepper sublist
190 pl->sublist("Default Integrator").set("Stepper Name", "General IMEX RK");
191 } else {
192 pl->sublist("Default Stepper").set("Stepper Type", stepperType);
193 }
194
195 // Set the step size
196 if (n == nTimeStepSizes-1) dt /= 10.0;
197 else dt /= 2;
198
199 // Setup the Integrator and reset initial time step
200 pl->sublist("Default Integrator")
201 .sublist("Time Step Control").set("Initial Time Step", dt);
202 pl->sublist("Default Integrator")
203 .sublist("Time Step Control").remove("Time Step Control Strategy");
204 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
205 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
206 order = integrator->getStepper()->getOrder();
207
208 // Integrate to timeMax
209 bool integratorStatus = integrator->advanceTime();
210 TEST_ASSERT(integratorStatus)
211
212 // Test if at 'Final Time'
213 double time = integrator->getTime();
214 double timeFinal =pl->sublist("Default Integrator")
215 .sublist("Time Step Control").get<double>("Final Time");
216 double tol = 100.0 * std::numeric_limits<double>::epsilon();
217 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
218
219 // Store off the final solution and step size
220 auto solution = Thyra::createMember(model->get_x_space());
221 auto sensitivity = Thyra::createMember(model->get_x_space());
222 Thyra::copy(*(integrator->getX()),solution.ptr());
223 Thyra::copy(*(integrator->getDxDp()->col(0)),sensitivity.ptr());
224 solutions.push_back(solution);
225 sensitivities.push_back(sensitivity);
226 StepSize.push_back(dt);
227
228 // Output finest temporal solution for plotting
229 if ((n == 0) || (n == nTimeStepSizes-1)) {
230 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
231
232 std::string fname = "Tempus_"+stepperName+"_VanDerPol_Sens-Ref.dat";
233 if (n == 0) fname = "Tempus_"+stepperName+"_VanDerPol_Sens.dat";
234 std::ofstream ftmp(fname);
235 RCP<const SolutionHistory<double> > solutionHistory =
236 integrator->getSolutionHistory();
237 int nStates = solutionHistory->getNumStates();
238 for (int i=0; i<nStates; i++) {
239 RCP<const SolutionState<double> > solutionState =
240 (*solutionHistory)[i];
241 RCP<const DMVPV> x_prod =
242 Teuchos::rcp_dynamic_cast<const DMVPV>(solutionState->getX());
243 RCP<const Thyra::VectorBase<double> > x =
244 x_prod->getMultiVector()->col(0);
245 RCP<const Thyra::VectorBase<double> > dxdp =
246 x_prod->getMultiVector()->col(1);
247 double ttime = solutionState->getTime();
248 ftmp << std::fixed << std::setprecision(7)
249 << ttime << " "
250 << std::setw(11) << get_ele(*x, 0) << " "
251 << std::setw(11) << get_ele(*x, 1) << " "
252 << std::setw(11) << get_ele(*dxdp, 0) << " "
253 << std::setw(11) << get_ele(*dxdp, 1)
254 << std::endl;
255 }
256 ftmp.close();
257 }
258 }
259
260 // Calculate the error - use the most temporally refined mesh for
261 // the reference solution.
262 auto ref_solution = solutions[solutions.size()-1];
263 auto ref_sensitivity = sensitivities[solutions.size()-1];
264 std::vector<double> StepSizeCheck;
265 for (std::size_t i=0; i < (solutions.size()-1); ++i) {
266 auto sol = solutions[i];
267 auto sen = sensitivities[i];
268 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
269 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
270 const double L2norm_sol = Thyra::norm_2(*sol);
271 const double L2norm_sen = Thyra::norm_2(*sen);
272 const double L2norm =
273 std::sqrt(L2norm_sol*L2norm_sol + L2norm_sen*L2norm_sen);
274 StepSizeCheck.push_back(StepSize[i]);
275 ErrorNorm.push_back(L2norm);
276
277 //*my_out << " n = " << i << " dt = " << StepSize[i]
278 // << " error = " << L2norm << std::endl;
279 }
280
281 // Check the order and intercept
282 double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
283 out << " Stepper = " << stepperType << std::endl;
284 out << " =========================" << std::endl;
285 out << " Expected order: " << order << std::endl;
286 out << " Observed order: " << slope << std::endl;
287 out << " =========================" << std::endl;
288 TEST_FLOATING_EQUALITY( slope, stepperOrders[m], 0.02 );
289 TEST_FLOATING_EQUALITY( ErrorNorm[0], stepperErrors[m], 1.0e-4 );
290
291 // Write error data
292 {
293 std::ofstream ftmp("Tempus_"+stepperName+"_VanDerPol_Sens-Error.dat");
294 double error0 = 0.8*ErrorNorm[0];
295 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
296 ftmp << StepSizeCheck[n] << " " << ErrorNorm[n] << " "
297 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
298 }
299 ftmp.close();
300 }
301 }
302 Teuchos::TimeMonitor::summarize();
303}
304
305
306} // namespace Tempus_Test
std::string method_name
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)