12 #include "Teuchos_DefaultComm.hpp"
14 #include "Tempus_config.hpp"
15 #include "Tempus_IntegratorBasic.hpp"
16 #include "Tempus_StepperTrapezoidal.hpp"
18 #include "../TestModels/SinCosModel.hpp"
19 #include "../TestModels/VanDerPolModel.hpp"
20 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
27 namespace Tempus_Test {
29 using Teuchos::getParametersFromXmlFile;
33 using Teuchos::rcp_const_cast;
34 using Teuchos::sublist;
46 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
57 Tempus::createIntegratorBasic<double>(tempusPL, model);
62 integrator->getStepper()->getValidParameters();
63 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
66 out <<
"stepperPL -------------- \n"
67 << *stepperPL << std::endl;
68 out <<
"defaultPL -------------- \n"
69 << *defaultPL << std::endl;
77 Tempus::createIntegratorBasic<double>(
78 model, std::string(
"Trapezoidal Method"));
82 integrator->getStepper()->getValidParameters();
84 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
87 out <<
"stepperPL -------------- \n"
88 << *stepperPL << std::endl;
89 out <<
"defaultPL -------------- \n"
90 << *defaultPL << std::endl;
101 std::vector<std::string> options;
102 options.push_back(
"Default Parameters");
103 options.push_back(
"ICConsistency and Check");
105 for (
const auto& option : options) {
108 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
131 stepper->setModel(model);
132 if (option ==
"ICConsistency and Check") {
133 stepper->setICConsistency(
"Consistent");
134 stepper->setICConsistencyCheck(
true);
136 stepper->initialize();
142 timeStepControl->setInitIndex(tscPL.
get<
int>(
"Initial Time Index"));
143 timeStepControl->setInitTime(tscPL.
get<
double>(
"Initial Time"));
144 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
145 timeStepControl->setInitTimeStep(dt);
146 timeStepControl->initialize();
149 auto inArgsIC = model->getNominalValues();
154 icState->setTime(timeStepControl->getInitTime());
155 icState->setIndex(timeStepControl->getInitIndex());
156 icState->setTimeStep(0.0);
157 icState->setOrder(stepper->getOrder());
162 solutionHistory->setName(
"Forward States");
164 solutionHistory->setStorageLimit(2);
165 solutionHistory->addState(icState);
169 Tempus::createIntegratorBasic<double>();
170 integrator->setStepper(stepper);
171 integrator->setTimeStepControl(timeStepControl);
172 integrator->setSolutionHistory(solutionHistory);
174 integrator->initialize();
177 bool integratorStatus = integrator->advanceTime();
181 double time = integrator->getTime();
182 double timeFinal = pl->
sublist(
"Default Integrator")
184 .
get<
double>(
"Final Time");
190 model->getExactSolution(time).get_x();
194 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
197 out <<
" Stepper = " << stepper->description() <<
"\n with "
198 << option << std::endl;
199 out <<
" =========================" << std::endl;
200 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
201 << get_ele(*(x_exact), 1) << std::endl;
202 out <<
" Computed solution: " << get_ele(*(x), 0) <<
" "
203 << get_ele(*(x), 1) << std::endl;
204 out <<
" Difference : " << get_ele(*(xdiff), 0) <<
" "
205 << get_ele(*(xdiff), 1) << std::endl;
206 out <<
" =========================" << std::endl;
207 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4);
208 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4);
217 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
218 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
219 std::vector<double> StepSize;
220 std::vector<double> xErrorNorm;
221 std::vector<double> xDotErrorNorm;
222 const int nTimeStepSizes = 7;
225 for (
int n = 0; n < nTimeStepSizes; n++) {
228 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
243 pl->
sublist(
"Default Integrator")
245 .
set(
"Initial Time Step", dt);
246 integrator = Tempus::createIntegratorBasic<double>(pl, model);
254 model->getNominalValues().get_x()->clone_v();
256 model->getNominalValues().get_x_dot()->clone_v();
257 integrator->initializeSolutionHistory(0.0, x0, xdot0);
261 bool integratorStatus = integrator->advanceTime();
265 time = integrator->getTime();
266 double timeFinal = pl->sublist(
"Default Integrator")
267 .sublist(
"Time Step Control")
268 .
get<
double>(
"Final Time");
274 integrator->getSolutionHistory();
275 writeSolution(
"Tempus_Trapezoidal_SinCos.dat", solutionHistory);
278 for (
int i = 0; i < solutionHistory->getNumStates(); i++) {
279 double time_i = (*solutionHistory)[i]->getTime();
282 model->getExactSolution(time_i).get_x()),
284 model->getExactSolution(time_i).get_x_dot()));
285 state->setTime((*solutionHistory)[i]->getTime());
286 solnHistExact->addState(state);
288 writeSolution(
"Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
292 StepSize.push_back(dt);
293 auto solution = Thyra::createMember(model->get_x_space());
294 Thyra::copy(*(integrator->getX()), solution.ptr());
295 solutions.push_back(solution);
296 auto solutionDot = Thyra::createMember(model->get_x_space());
297 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
298 solutionsDot.push_back(solutionDot);
299 if (n == nTimeStepSizes - 1) {
300 StepSize.push_back(0.0);
301 auto solutionExact = Thyra::createMember(model->get_x_space());
302 Thyra::copy(*(model->getExactSolution(time).get_x()),
303 solutionExact.ptr());
304 solutions.push_back(solutionExact);
305 auto solutionDotExact = Thyra::createMember(model->get_x_space());
306 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
307 solutionDotExact.ptr());
308 solutionsDot.push_back(solutionDotExact);
313 double xDotSlope = 0.0;
315 double order = stepper->getOrder();
316 writeOrderError(
"Tempus_Trapezoidal_SinCos-Error.dat", stepper, StepSize,
317 solutions, xErrorNorm, xSlope, solutionsDot, xDotErrorNorm,
333 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
334 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
335 std::vector<double> StepSize;
336 std::vector<double> xErrorNorm;
337 std::vector<double> xDotErrorNorm;
338 const int nTimeStepSizes = 4;
341 for (
int n = 0; n < nTimeStepSizes; n++) {
344 getParametersFromXmlFile(
"Tempus_Trapezoidal_VanDerPol.xml");
352 if (n == nTimeStepSizes - 1) dt /= 10.0;
358 .
set(
"Initial Time Step", dt);
359 integrator = Tempus::createIntegratorBasic<double>(pl, model);
362 bool integratorStatus = integrator->advanceTime();
366 time = integrator->getTime();
367 double timeFinal = pl->sublist(
"Demo Integrator")
368 .sublist(
"Time Step Control")
369 .
get<
double>(
"Final Time");
370 double tol = 100.0 * std::numeric_limits<double>::epsilon();
374 StepSize.push_back(dt);
375 auto solution = Thyra::createMember(model->get_x_space());
376 Thyra::copy(*(integrator->getX()), solution.ptr());
377 solutions.push_back(solution);
378 auto solutionDot = Thyra::createMember(model->get_x_space());
379 Thyra::copy(*(integrator->getXDot()), solutionDot.ptr());
380 solutionsDot.push_back(solutionDot);
384 if ((n == 0) || (n == nTimeStepSizes - 1)) {
385 std::string fname =
"Tempus_Trapezoidal_VanDerPol-Ref.dat";
386 if (n == 0) fname =
"Tempus_Trapezoidal_VanDerPol.dat";
388 integrator->getSolutionHistory();
394 double xDotSlope = 0.0;
396 double order = stepper->getOrder();
397 writeOrderError(
"Tempus_Trapezoidal_VanDerPol-Error.dat", stepper, StepSize,
398 solutions, xErrorNorm, xSlope, solutionsDot, xDotErrorNorm,
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
Trapezoidal method time stepper.
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)
void writeOrderError(const std::string filename, Teuchos::RCP< Tempus::Stepper< Scalar >> stepper, std::vector< Scalar > &StepSize, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutions, std::vector< Scalar > &xErrorNorm, Scalar &xSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDot, std::vector< Scalar > &xDotErrorNorm, Scalar &xDotSlope, std::vector< Teuchos::RCP< Thyra::VectorBase< Scalar >>> &solutionsDotDot, std::vector< Scalar > &xDotDotErrorNorm, Scalar &xDotDotSlope, Teuchos::FancyOStream &out)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar >> solutionHistory)
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)
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Keep a fix number of states.
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.