12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperBDF2.hpp"
18 #include "../TestModels/VanDerPolModel.hpp"
19 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
26 namespace Tempus_Test {
28 using Teuchos::getParametersFromXmlFile;
32 using Teuchos::rcp_const_cast;
33 using Teuchos::sublist;
43 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
44 std::vector<double> StepSize;
45 std::vector<double> ErrorNorm;
49 getParametersFromXmlFile(
"Tempus_BDF2_VanDerPol.xml");
54 double dt = pl->
sublist(
"Demo Integrator")
56 .
get<
double>(
"Initial Time Step");
60 const int nTimeStepSizes = vdpm_pl->
get<
int>(
"Number of Time Step Sizes", 3);
64 for (
int n = 0; n < nTimeStepSizes; n++) {
70 if (n == nTimeStepSizes - 1) dt /= 10.0;
75 .
set(
"Initial Time Step", dt);
77 Tempus::createIntegratorBasic<double>(pl, model);
78 order = integrator->getStepper()->getOrder();
81 bool integratorStatus = integrator->advanceTime();
85 double time = integrator->getTime();
86 double timeFinal = pl->sublist(
"Demo Integrator")
87 .sublist(
"Time Step Control")
88 .
get<
double>(
"Final Time");
89 double tol = 100.0 * std::numeric_limits<double>::epsilon();
93 auto solution = Thyra::createMember(model->get_x_space());
94 Thyra::copy(*(integrator->getX()), solution.ptr());
95 solutions.push_back(solution);
96 StepSize.push_back(dt);
100 if ((n == 0) || (n == nTimeStepSizes - 1)) {
101 std::string fname =
"Tempus_BDF2_VanDerPol-Ref.dat";
102 if (n == 0) fname =
"Tempus_BDF2_VanDerPol.dat";
103 std::ofstream ftmp(fname);
105 integrator->getSolutionHistory();
106 int nStates = solutionHistory->getNumStates();
107 for (
int i = 0; i < nStates; i++) {
110 double ttime = solutionState->getTime();
111 ftmp << ttime <<
" " << get_ele(*x, 0) <<
" " << get_ele(*x, 1)
120 auto ref_solution = solutions[solutions.size() - 1];
121 std::vector<double> StepSizeCheck;
122 for (std::size_t i = 0; i < (solutions.size() - 1); ++i) {
123 auto tmp = solutions[i];
124 Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
125 const double L2norm = Thyra::norm_2(*tmp);
126 StepSizeCheck.push_back(StepSize[i]);
127 ErrorNorm.push_back(L2norm);
130 if (nTimeStepSizes > 2) {
133 computeLinearRegressionLogLog<double>(StepSizeCheck, ErrorNorm);
134 out <<
" Stepper = BDF2" << std::endl;
135 out <<
" =========================" << std::endl;
136 out <<
" Expected order: " << order << std::endl;
137 out <<
" Observed order: " << slope << std::endl;
138 out <<
" =========================" << std::endl;
140 out <<
"\n\n ** Slope on BDF2 Method = " << slope <<
"\n"
146 std::ofstream ftmp(
"Tempus_BDF2_VanDerPol-Error.dat");
147 double error0 = 0.8 * ErrorNorm[0];
148 for (std::size_t n = 0; n < StepSizeCheck.size(); n++) {
149 ftmp << StepSizeCheck[n] <<
" " << ErrorNorm[n] <<
" "
150 << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)
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...
van der Pol model problem for nonlinear electrical circuit.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Solution state for integrators and steppers.