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"
31 namespace Tempus_Test {
33 using Teuchos::getParametersFromXmlFile;
36 using Teuchos::sublist;
48 std::vector<std::string> stepperTypes;
49 stepperTypes.push_back(
"IMEX RK 1st order");
50 stepperTypes.push_back(
"IMEX RK SSP2");
51 stepperTypes.push_back(
"IMEX RK ARS 233");
52 stepperTypes.push_back(
"General IMEX RK");
54 std::vector<double> stepperOrders;
55 std::vector<double> stepperErrors;
56 if (use_combined_method) {
57 stepperOrders.push_back(1.1198);
58 stepperOrders.push_back(1.98931);
59 stepperOrders.push_back(2.60509);
60 stepperOrders.push_back(1.992);
62 stepperErrors.push_back(0.00619674);
63 stepperErrors.push_back(0.294989);
64 stepperErrors.push_back(0.0062125);
65 stepperErrors.push_back(0.142489);
68 stepperOrders.push_back(1.1198);
69 stepperOrders.push_back(2.05232);
70 stepperOrders.push_back(2.43013);
71 stepperOrders.push_back(2.07975);
73 stepperErrors.push_back(0.00619674);
74 stepperErrors.push_back(0.0534172);
75 stepperErrors.push_back(0.00224533);
76 stepperErrors.push_back(0.032632);
78 std::vector<double> stepperInitDt;
79 stepperInitDt.push_back(0.0125);
80 stepperInitDt.push_back(0.05);
81 stepperInitDt.push_back(0.05);
82 stepperInitDt.push_back(0.05);
87 std::vector<std::string>::size_type m;
88 for (m = 0; m != stepperTypes.size(); m++) {
89 std::string stepperType = stepperTypes[m];
90 std::string stepperName = stepperTypes[m];
91 std::replace(stepperName.begin(), stepperName.end(),
' ',
'_');
92 std::replace(stepperName.begin(), stepperName.end(),
'/',
'.');
94 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
95 std::vector<RCP<Thyra::VectorBase<double>>> sensitivities;
96 std::vector<double> StepSize;
97 std::vector<double> ErrorNorm;
98 const int nTimeStepSizes = 3;
99 double dt = stepperInitDt[m];
101 for (
int n = 0; n < nTimeStepSizes; n++) {
104 getParametersFromXmlFile(
"Tempus_IMEX_RK_VanDerPol.xml");
108 vdpmPL->
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
119 explicitModel, implicitModel));
124 if (use_combined_method)
125 sens_pl.
set(
"Sensitivity Method",
"Combined");
127 sens_pl.
set(
"Sensitivity Method",
"Staggered");
128 sens_pl.
set(
"Reuse State Linear Solver",
true);
130 sens_pl.
set(
"Use DfDp as Tangent", use_dfdp_as_tangent);
134 interp_pl.
set(
"Interpolator Type",
"Lagrange");
135 interp_pl.
set(
"Order", 2);
138 if (stepperType ==
"General IMEX RK") {
140 pl->
sublist(
"Default Integrator")
141 .
set(
"Stepper Name",
"General IMEX RK");
144 pl->
sublist(
"Default Stepper").
set(
"Stepper Type", stepperType);
148 if (n == nTimeStepSizes - 1)
154 pl->
sublist(
"Default Integrator")
156 .
set(
"Initial Time Step", dt);
157 pl->
sublist(
"Default Integrator")
159 .
remove(
"Time Step Control Strategy");
161 Tempus::createIntegratorForwardSensitivity<double>(pl, model);
162 order = integrator->getStepper()->getOrder();
165 bool integratorStatus = integrator->advanceTime();
169 double time = integrator->getTime();
170 double timeFinal = pl->sublist(
"Default Integrator")
171 .sublist(
"Time Step Control")
172 .
get<
double>(
"Final Time");
173 double tol = 100.0 * std::numeric_limits<double>::epsilon();
177 auto solution = Thyra::createMember(model->get_x_space());
178 auto sensitivity = Thyra::createMember(model->get_x_space());
179 Thyra::copy(*(integrator->getX()), solution.ptr());
180 Thyra::copy(*(integrator->getDxDp()->col(0)), sensitivity.ptr());
181 solutions.push_back(solution);
182 sensitivities.push_back(sensitivity);
183 StepSize.push_back(dt);
186 if (comm->getRank() == 0 && ((n == 0) || (n == nTimeStepSizes - 1))) {
189 std::string fname =
"Tempus_" + stepperName +
"_VanDerPol_Sens-Ref.dat";
190 if (n == 0) fname =
"Tempus_" + stepperName +
"_VanDerPol_Sens.dat";
191 std::ofstream ftmp(fname);
193 integrator->getSolutionHistory();
194 int nStates = solutionHistory->getNumStates();
195 for (
int i = 0; i < nStates; i++) {
197 (*solutionHistory)[i];
199 Teuchos::rcp_dynamic_cast<
const DMVPV>(solutionState->getX());
201 x_prod->getMultiVector()->col(0);
203 x_prod->getMultiVector()->col(1);
204 double ttime = solutionState->getTime();
205 ftmp << std::fixed << std::setprecision(7) << ttime <<
" "
206 << std::setw(11) << get_ele(*x, 0) <<
" " << std::setw(11)
207 << get_ele(*x, 1) <<
" " << std::setw(11) << get_ele(*dxdp, 0)
208 <<
" " << std::setw(11) << get_ele(*dxdp, 1) << std::endl;
216 auto ref_solution = solutions[solutions.size() - 1];
217 auto ref_sensitivity = sensitivities[solutions.size() - 1];
218 std::vector<double> StepSizeCheck;
219 for (std::size_t i = 0; i < (solutions.size() - 1); ++i) {
220 auto sol = solutions[i];
221 auto sen = sensitivities[i];
222 Thyra::Vp_StV(sol.ptr(), -1.0, *ref_solution);
223 Thyra::Vp_StV(sen.ptr(), -1.0, *ref_sensitivity);
224 const double L2norm_sol = Thyra::norm_2(*sol);
225 const double L2norm_sen = Thyra::norm_2(*sen);
226 const double L2norm =
227 std::sqrt(L2norm_sol * L2norm_sol + L2norm_sen * L2norm_sen);
228 StepSizeCheck.push_back(StepSize[i]);
229 ErrorNorm.push_back(L2norm);
231 out <<
" n = " << i <<
" dt = " << StepSize[i] <<
" error = " << L2norm
237 computeLinearRegressionLogLog<double>(StepSizeCheck, ErrorNorm);
238 out <<
" Stepper = " << stepperType << std::endl;
239 out <<
" =========================" << std::endl;
240 out <<
" Expected order: " << order << std::endl;
241 out <<
" Observed order: " << slope << std::endl;
242 out <<
" =========================" << std::endl;
248 std::ofstream ftmp(
"Tempus_" + stepperName +
"_VanDerPol_Sens-Error.dat");
249 double error0 = 0.8 * ErrorNorm[0];
250 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
251 ftmp << StepSizeCheck[n] <<
" " << ErrorNorm[n] <<
" "
252 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
van der Pol model formulated for IMEX-RK.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
ModelEvaluator pair for implicit and explicit (IMEX) evaulations.
#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)
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.