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"
25#include "../TestUtils/Tempus_ConvergenceTestUtils.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.
van der Pol model formulated for the partitioned IMEX-RK.
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)