9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12#include "Teuchos_DefaultComm.hpp"
14#include "Thyra_VectorStdOps.hpp"
15#include "Thyra_MultiVectorStdOps.hpp"
16#include "Thyra_DefaultMultiVectorProductVector.hpp"
17#include "Thyra_DefaultProductVector.hpp"
19#include "Tempus_IntegratorBasic.hpp"
20#include "Tempus_IntegratorForwardSensitivity.hpp"
21#include "Tempus_WrapperModelEvaluatorPairPartIMEX_Basic.hpp"
23#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
24#include "../TestModels/VanDerPol_IMEXPart_ImplicitModel.hpp"
25#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
33using Teuchos::ParameterList;
34using Teuchos::sublist;
35using Teuchos::getParametersFromXmlFile;
45 const bool use_combined_method,
46 const bool use_dfdp_as_tangent,
47 Teuchos::FancyOStream &out,
bool &success)
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" );
57 auto it = std::find(stepperTypes.begin(), stepperTypes.end(),
method_name);
58 TEUCHOS_TEST_FOR_EXCEPTION(it == stepperTypes.end(), std::logic_error,
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);
71 stepperErrors.push_back(0.00820931);
72 stepperErrors.push_back(0.287112);
73 stepperErrors.push_back(0.00646096);
74 stepperErrors.push_back(0.148848);
77 stepperOrders.push_back(1.07932);
78 stepperOrders.push_back(1.97396);
79 stepperOrders.push_back(2.63724);
80 stepperOrders.push_back(1.99133);
82 stepperErrors.push_back(0.055626);
83 stepperErrors.push_back(0.198898);
84 stepperErrors.push_back(0.00614135);
85 stepperErrors.push_back(0.0999881);
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);
95 stepperErrors.push_back(0.00619674);
96 stepperErrors.push_back(0.294989);
97 stepperErrors.push_back(0.0062125);
98 stepperErrors.push_back(0.142489);
101 stepperOrders.push_back(1.07932);
102 stepperOrders.push_back(1.97396);
103 stepperOrders.push_back(2.63724);
104 stepperOrders.push_back(1.99133);
106 stepperErrors.push_back(0.055626);
107 stepperErrors.push_back(0.198898);
108 stepperErrors.push_back(0.00614135);
109 stepperErrors.push_back(0.0999881);
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);
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);
126 std::vector<std::string>::size_type m;
127 for(m = 0; m != stepperTypes.size(); m++) {
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(),
'/',
'.');
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;
143 double dt = stepperInitDt[m];
145 for (
int n=0; n<nTimeStepSizes; n++) {
148 RCP<ParameterList> pList =
149 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
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 =
160 RCP<VanDerPol_IMEXPart_ImplicitModel<double> > implicitModel =
164 const int numExplicitBlocks = 1;
165 const int parameterIndex = 4;
166 RCP<Tempus::WrapperModelEvaluatorPairPartIMEX_Basic<double> > model =
169 explicitModel, implicitModel,
170 numExplicitBlocks, parameterIndex));
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");
178 sens_pl.set(
"Sensitivity Method",
"Staggered");
179 sens_pl.set(
"Reuse State Linear Solver",
true);
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);
188 if (stepperType ==
"General Partitioned IMEX RK"){
190 pl->sublist(
"Default Integrator").set(
"Stepper Name",
"General IMEX RK");
192 pl->sublist(
"Default Stepper").set(
"Stepper Type", stepperType);
196 if (n == nTimeStepSizes-1) dt /= 10.0;
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();
209 bool integratorStatus = integrator->advanceTime();
210 TEST_ASSERT(integratorStatus)
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);
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);
229 if ((n == 0) || (n == nTimeStepSizes-1)) {
230 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
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)
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)
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);
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 );
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;
302 Teuchos::TimeMonitor::summarize();
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.
van der Pol model formulated for IMEX.
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)