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"
30 namespace Tempus_Test {
33 using Teuchos::ParameterList;
34 using Teuchos::sublist;
35 using 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" );
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,
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++) {
130 if (method_name !=
"" && stepperTypes[m] != method_name)
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 RCP<Tempus::IntegratorForwardSensitivity<double> > integrator =
203 Tempus::integratorForwardSensitivity<double>(pl, model);
204 order = integrator->getStepper()->getOrder();
207 bool integratorStatus = integrator->advanceTime();
208 TEST_ASSERT(integratorStatus)
211 double time = integrator->getTime();
212 double timeFinal =pl->sublist(
"Default Integrator")
213 .sublist(
"Time Step Control").get<
double>(
"Final Time");
214 double tol = 100.0 * std::numeric_limits<double>::epsilon();
215 TEST_FLOATING_EQUALITY(time, timeFinal, tol);
218 auto solution = Thyra::createMember(model->get_x_space());
219 auto sensitivity = Thyra::createMember(model->get_x_space());
220 Thyra::copy(*(integrator->getX()),solution.ptr());
221 Thyra::copy(*(integrator->getDxDp()->col(0)),sensitivity.ptr());
222 solutions.push_back(solution);
223 sensitivities.push_back(sensitivity);
224 StepSize.push_back(dt);
227 if ((n == 0) or (n == nTimeStepSizes-1)) {
228 typedef Thyra::DefaultMultiVectorProductVector<double> DMVPV;
230 std::string fname =
"Tempus_"+stepperName+
"_VanDerPol_Sens-Ref.dat";
231 if (n == 0) fname =
"Tempus_"+stepperName+
"_VanDerPol_Sens.dat";
232 std::ofstream ftmp(fname);
234 integrator->getSolutionHistory();
235 int nStates = solutionHistory->getNumStates();
236 for (
int i=0; i<nStates; i++) {
237 RCP<const SolutionState<double> > solutionState =
238 (*solutionHistory)[i];
239 RCP<const DMVPV> x_prod =
240 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
241 RCP<const Thyra::VectorBase<double> > x =
242 x_prod->getMultiVector()->col(0);
243 RCP<const Thyra::VectorBase<double> > dxdp =
244 x_prod->getMultiVector()->col(1);
245 double ttime = solutionState->getTime();
246 ftmp << std::fixed << std::setprecision(7)
248 << std::setw(11) << get_ele(*x, 0) <<
" "
249 << std::setw(11) << get_ele(*x, 1) <<
" "
250 << std::setw(11) << get_ele(*dxdp, 0) <<
" "
251 << std::setw(11) << get_ele(*dxdp, 1)
260 auto ref_solution = solutions[solutions.size()-1];
261 auto ref_sensitivity = sensitivities[solutions.size()-1];
262 std::vector<double> StepSizeCheck;
263 for (std::size_t i=0; i < (solutions.size()-1); ++i) {
264 auto sol = solutions[i];
265 auto sen = sensitivities[i];
266 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
267 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
268 const double L2norm_sol = Thyra::norm_2(*sol);
269 const double L2norm_sen = Thyra::norm_2(*sen);
270 const double L2norm =
271 std::sqrt(L2norm_sol*L2norm_sol + L2norm_sen*L2norm_sen);
272 StepSizeCheck.push_back(StepSize[i]);
273 ErrorNorm.push_back(L2norm);
275 *my_out <<
" n = " << i <<
" dt = " << StepSize[i]
276 <<
" error = " << L2norm << std::endl;
280 double slope = computeLinearRegressionLogLog<double>(StepSizeCheck,ErrorNorm);
281 std::cout <<
" Stepper = " << stepperType << std::endl;
282 std::cout <<
" =========================" << std::endl;
283 std::cout <<
" Expected order: " << order << std::endl;
284 std::cout <<
" Observed order: " << slope << std::endl;
285 std::cout <<
" =========================" << std::endl;
286 TEST_FLOATING_EQUALITY( slope, stepperOrders[m], 0.02 );
287 TEST_FLOATING_EQUALITY( ErrorNorm[0], stepperErrors[m], 1.0e-4 );
291 std::ofstream ftmp(
"Tempus_"+stepperName+
"_VanDerPol_Sens-Error.dat");
292 double error0 = 0.8*ErrorNorm[0];
293 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
294 ftmp << StepSizeCheck[n] <<
" " << ErrorNorm[n] <<
" "
295 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
300 Teuchos::TimeMonitor::summarize();
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)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
van der Pol model formulated for IMEX.
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...