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 {
32 using Teuchos::getParametersFromXmlFile;
35 using Teuchos::sublist;
44 const bool use_combined_method,
48 std::vector<std::string> stepperTypes;
49 stepperTypes.push_back(
"Partitioned IMEX RK 1st order");
50 stepperTypes.push_back(
"Partitioned IMEX RK SSP2");
51 stepperTypes.push_back(
"Partitioned IMEX RK ARS 233");
52 stepperTypes.push_back(
"General Partitioned IMEX RK");
55 if (method_name !=
"") {
56 auto it = std::find(stepperTypes.begin(), stepperTypes.end(),
method_name);
61 std::vector<double> stepperOrders;
62 std::vector<double> stepperErrors;
63 if (use_dfdp_as_tangent) {
64 if (use_combined_method) {
65 stepperOrders.push_back(1.16082);
66 stepperOrders.push_back(1.97231);
67 stepperOrders.push_back(2.5914);
68 stepperOrders.push_back(1.99148);
70 stepperErrors.push_back(0.00820931);
71 stepperErrors.push_back(0.287112);
72 stepperErrors.push_back(0.00646096);
73 stepperErrors.push_back(0.148848);
76 stepperOrders.push_back(1.07932);
77 stepperOrders.push_back(1.97396);
78 stepperOrders.push_back(2.63724);
79 stepperOrders.push_back(1.99133);
81 stepperErrors.push_back(0.055626);
82 stepperErrors.push_back(0.198898);
83 stepperErrors.push_back(0.00614135);
84 stepperErrors.push_back(0.0999881);
88 if (use_combined_method) {
89 stepperOrders.push_back(1.1198);
90 stepperOrders.push_back(1.98931);
91 stepperOrders.push_back(2.60509);
92 stepperOrders.push_back(1.992);
94 stepperErrors.push_back(0.00619674);
95 stepperErrors.push_back(0.294989);
96 stepperErrors.push_back(0.0062125);
97 stepperErrors.push_back(0.142489);
100 stepperOrders.push_back(1.07932);
101 stepperOrders.push_back(1.97396);
102 stepperOrders.push_back(2.63724);
103 stepperOrders.push_back(1.99133);
105 stepperErrors.push_back(0.055626);
106 stepperErrors.push_back(0.198898);
107 stepperErrors.push_back(0.00614135);
108 stepperErrors.push_back(0.0999881);
112 std::vector<double> stepperInitDt;
113 stepperInitDt.push_back(0.0125);
114 stepperInitDt.push_back(0.05);
115 stepperInitDt.push_back(0.05);
116 stepperInitDt.push_back(0.05);
121 std::vector<std::string>::size_type m;
122 for (m = 0; m != stepperTypes.size(); m++) {
124 if (method_name !=
"" && stepperTypes[m] != method_name)
continue;
126 std::string stepperType = stepperTypes[m];
127 std::string stepperName = stepperTypes[m];
128 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
129 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
131 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
132 std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
133 std::vector<double> StepSize;
134 std::vector<double> ErrorNorm;
135 const int nTimeStepSizes = 3;
136 double dt = stepperInitDt[m];
138 for (
int n = 0; n < nTimeStepSizes; n++) {
141 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
145 vdpmPL->
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
146 const bool useProductVector =
true;
155 const int numExplicitBlocks = 1;
156 const int parameterIndex = 4;
160 explicitModel, implicitModel, numExplicitBlocks,
166 if (use_combined_method)
167 sens_pl.
set(
"Sensitivity Method",
"Combined");
169 sens_pl.
set(
"Sensitivity Method",
"Staggered");
170 sens_pl.
set(
"Reuse State Linear Solver",
true);
172 sens_pl.
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
176 interp_pl.
set(
"Interpolator Type",
"Lagrange");
177 interp_pl.
set(
"Order", 2);
180 if (stepperType ==
"General Partitioned IMEX RK") {
182 pl->
sublist(
"Default Integrator")
183 .
set(
"Stepper Name",
"General IMEX RK");
186 pl->
sublist(
"Default Stepper").
set(
"Stepper Type", stepperType);
190 if (n == nTimeStepSizes - 1)
196 pl->
sublist(
"Default Integrator")
198 .
set(
"Initial Time Step", dt);
199 pl->
sublist(
"Default Integrator")
201 .
remove(
"Time Step Control Strategy");
203 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
204 order = integrator->getStepper()->getOrder();
207 bool integratorStatus = integrator->advanceTime();
211 double time = integrator->getTime();
212 double timeFinal = pl->sublist(
"Default Integrator")
213 .sublist(
"Time Step Control")
214 .
get<
double>(
"Final Time");
215 double tol = 100.0 * std::numeric_limits<double>::epsilon();
219 auto solution = Thyra::createMember(model->get_x_space());
220 auto sensitivity = Thyra::createMember(model->get_x_space());
221 Thyra::copy(*(integrator->getX()), solution.ptr());
222 Thyra::copy(*(integrator->getDxDp()->col(0)), sensitivity.ptr());
223 solutions.push_back(solution);
224 sensitivities.push_back(sensitivity);
225 StepSize.push_back(dt);
228 if ((n == 0) || (n == nTimeStepSizes - 1)) {
231 std::string fname =
"Tempus_" + stepperName +
"_VanDerPol_Sens-Ref.dat";
232 if (n == 0) fname =
"Tempus_" + stepperName +
"_VanDerPol_Sens.dat";
233 std::ofstream ftmp(fname);
235 integrator->getSolutionHistory();
236 int nStates = solutionHistory->getNumStates();
237 for (
int i = 0; i < nStates; i++) {
239 (*solutionHistory)[i];
241 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
243 x_prod->getMultiVector()->col(0);
245 x_prod->getMultiVector()->col(1);
246 double ttime = solutionState->getTime();
247 ftmp << std::fixed << std::setprecision(7) << ttime <<
" "
248 << std::setw(11) << get_ele(*x, 0) <<
" " << std::setw(11)
249 << get_ele(*x, 1) <<
" " << std::setw(11) << get_ele(*dxdp, 0)
250 <<
" " << std::setw(11) << get_ele(*dxdp, 1) << std::endl;
258 auto ref_solution = solutions[solutions.size() - 1];
259 auto ref_sensitivity = sensitivities[solutions.size() - 1];
260 std::vector<double> StepSizeCheck;
261 for (std::size_t i = 0; i < (solutions.size() - 1); ++i) {
262 auto sol = solutions[i];
263 auto sen = sensitivities[i];
264 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
265 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
266 const double L2norm_sol = Thyra::norm_2(*sol);
267 const double L2norm_sen = Thyra::norm_2(*sen);
268 const double L2norm =
269 std::sqrt(L2norm_sol * L2norm_sol + L2norm_sen * L2norm_sen);
270 StepSizeCheck.push_back(StepSize[i]);
271 ErrorNorm.push_back(L2norm);
279 computeLinearRegressionLogLog<double>(StepSizeCheck, ErrorNorm);
280 out <<
" Stepper = " << stepperType << std::endl;
281 out <<
" =========================" << std::endl;
282 out <<
" Expected order: " << order << std::endl;
283 out <<
" Observed order: " << slope << std::endl;
284 out <<
" =========================" << std::endl;
290 std::ofstream ftmp(
"Tempus_" + stepperName +
"_VanDerPol_Sens-Error.dat");
291 double error0 = 0.8 * ErrorNorm[0];
292 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
293 ftmp << StepSizeCheck[n] <<
" " << ErrorNorm[n] <<
" "
294 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
van der Pol model formulated for the partitioned IMEX-RK.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
void test_vdp_fsa(const bool use_combined_method, const bool use_dfdp_as_tangent, Teuchos::FancyOStream &out, bool &success)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
van der Pol model formulated for IMEX.
Solution state for integrators and steppers.