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_WrapperModelEvaluatorPairIMEX_Basic.hpp"
23#include "../TestModels/VanDerPolModel.hpp"
24#include "../TestModels/VanDerPol_IMEX_ExplicitModel.hpp"
25#include "../TestModels/VanDerPol_IMEX_ImplicitModel.hpp"
26#include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
34using Teuchos::ParameterList;
35using Teuchos::sublist;
36using Teuchos::getParametersFromXmlFile;
46 const bool use_dfdp_as_tangent,
47 Teuchos::FancyOStream &out,
bool &success)
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" );
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);
63 stepperErrors.push_back(0.00619674);
64 stepperErrors.push_back(0.294989);
65 stepperErrors.push_back(0.0062125);
66 stepperErrors.push_back(0.142489);
69 stepperOrders.push_back(1.1198);
70 stepperOrders.push_back(2.05232);
71 stepperOrders.push_back(2.43013);
72 stepperOrders.push_back(2.07975);
74 stepperErrors.push_back(0.00619674);
75 stepperErrors.push_back(0.0534172);
76 stepperErrors.push_back(0.00224533);
77 stepperErrors.push_back(0.032632);
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);
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);
92 std::vector<std::string>::size_type m;
93 for(m = 0; m != stepperTypes.size(); m++) {
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(),
'/',
'.');
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;
105 double dt = stepperInitDt[m];
107 for (
int n=0; n<nTimeStepSizes; n++) {
110 RCP<ParameterList> pList =
111 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
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 =
120 RCP<VanDerPol_IMEX_ImplicitModel<double> > implicitModel =
124 RCP<Tempus::WrapperModelEvaluatorPairIMEX_Basic<double> > model =
126 explicitModel, implicitModel));
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");
134 sens_pl.set(
"Sensitivity Method",
"Staggered");
135 sens_pl.set(
"Reuse State Linear Solver",
true);
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);
144 if (stepperType ==
"General IMEX RK"){
146 pl->sublist(
"Default Integrator").set(
"Stepper Name",
"General IMEX RK");
148 pl->sublist(
"Default Stepper").set(
"Stepper Type", stepperType);
152 if (n == nTimeStepSizes-1) dt /= 10.0;
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();
165 bool integratorStatus = integrator->advanceTime();
166 TEST_ASSERT(integratorStatus)
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);
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);
185 if (comm->getRank() == 0 && ((n == 0) || (n == nTimeStepSizes-1))) {
186 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
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)
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)
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);
233 *my_out <<
" n = " << i <<
" dt = " << StepSize[i]
234 <<
" error = " << L2norm << std::endl;
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 );
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;
258 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 IMEX.
van der Pol model formulated for IMEX-RK.
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)