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"
22 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
23 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
31 namespace Tempus_Test {
35 using Teuchos::rcp_const_cast;
37 using Teuchos::sublist;
38 using Teuchos::getParametersFromXmlFile;
51 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
62 Tempus::createIntegratorBasic<double>(tempusPL, model);
67 integrator->getStepper()->getValidParameters();
68 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
71 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
72 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
80 Tempus::createIntegratorBasic<double>(model, std::string(
"Trapezoidal Method"));
84 integrator->getStepper()->getValidParameters();
86 bool pass = haveSameValuesSorted(*stepperPL, *defaultPL,
true);
89 out <<
"stepperPL -------------- \n" << *stepperPL << std::endl;
90 out <<
"defaultPL -------------- \n" << *defaultPL << std::endl;
102 std::vector<std::string> options;
103 options.push_back(
"Default Parameters");
104 options.push_back(
"ICConsistency and Check");
106 for(
const auto& option: options) {
110 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
133 stepper->setModel(model);
134 if ( option ==
"ICConsistency and Check") {
135 stepper->setICConsistency(
"Consistent");
136 stepper->setICConsistencyCheck(
true);
138 stepper->initialize();
144 timeStepControl->setInitIndex(tscPL.
get<
int> (
"Initial Time Index"));
145 timeStepControl->setInitTime (tscPL.
get<
double>(
"Initial Time"));
146 timeStepControl->setFinalTime(tscPL.
get<
double>(
"Final Time"));
147 timeStepControl->setInitTimeStep(dt);
148 timeStepControl->initialize();
151 auto inArgsIC = model->getNominalValues();
156 icState->setTime (timeStepControl->getInitTime());
157 icState->setIndex (timeStepControl->getInitIndex());
158 icState->setTimeStep(0.0);
159 icState->setOrder (stepper->getOrder());
164 solutionHistory->setName(
"Forward States");
166 solutionHistory->setStorageLimit(2);
167 solutionHistory->addState(icState);
171 Tempus::createIntegratorBasic<double>();
172 integrator->setStepper(stepper);
173 integrator->setTimeStepControl(timeStepControl);
174 integrator->setSolutionHistory(solutionHistory);
176 integrator->initialize();
180 bool integratorStatus = integrator->advanceTime();
185 double time = integrator->getTime();
186 double timeFinal =pl->
sublist(
"Default Integrator")
187 .
sublist(
"Time Step Control").
get<
double>(
"Final Time");
193 model->getExactSolution(time).get_x();
197 Thyra::V_StVpStV(xdiff.
ptr(), 1.0, *x_exact, -1.0, *(x));
200 out <<
" Stepper = " << stepper->description()
201 <<
"\n with " << option << std::endl;
202 out <<
" =========================" << std::endl;
203 out <<
" Exact solution : " << get_ele(*(x_exact), 0) <<
" "
204 << get_ele(*(x_exact), 1) << std::endl;
205 out <<
" Computed solution: " << get_ele(*(x ), 0) <<
" "
206 << get_ele(*(x ), 1) << std::endl;
207 out <<
" Difference : " << get_ele(*(xdiff ), 0) <<
" "
208 << get_ele(*(xdiff ), 1) << std::endl;
209 out <<
" =========================" << std::endl;
210 TEST_FLOATING_EQUALITY(get_ele(*(x), 0), 0.841021, 1.0e-4 );
211 TEST_FLOATING_EQUALITY(get_ele(*(x), 1), 0.541002, 1.0e-4 );
221 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
222 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
223 std::vector<double> StepSize;
224 std::vector<double> xErrorNorm;
225 std::vector<double> xDotErrorNorm;
226 const int nTimeStepSizes = 7;
229 for (
int n=0; n<nTimeStepSizes; n++) {
233 getParametersFromXmlFile(
"Tempus_Trapezoidal_SinCos.xml");
248 pl->
sublist(
"Default Integrator")
249 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
250 integrator = Tempus::createIntegratorBasic<double>(pl, model);
258 model->getNominalValues().get_x()->clone_v();
260 model->getNominalValues().get_x_dot()->clone_v();
261 integrator->initializeSolutionHistory(0.0, x0, xdot0);
265 bool integratorStatus = integrator->advanceTime();
269 time = integrator->getTime();
270 double timeFinal =pl->sublist(
"Default Integrator")
271 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
277 integrator->getSolutionHistory();
278 writeSolution(
"Tempus_Trapezoidal_SinCos.dat", solutionHistory);
281 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
282 double time_i = (*solutionHistory)[i]->getTime();
285 model->getExactSolution(time_i).get_x()),
287 model->getExactSolution(time_i).get_x_dot()));
288 state->setTime((*solutionHistory)[i]->getTime());
289 solnHistExact->addState(state);
291 writeSolution(
"Tempus_Trapezoidal_SinCos-Ref.dat", solnHistExact);
295 StepSize.push_back(dt);
296 auto solution = Thyra::createMember(model->get_x_space());
297 Thyra::copy(*(integrator->getX()),solution.ptr());
298 solutions.push_back(solution);
299 auto solutionDot = Thyra::createMember(model->get_x_space());
300 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
301 solutionsDot.push_back(solutionDot);
302 if (n == nTimeStepSizes-1) {
303 StepSize.push_back(0.0);
304 auto solutionExact = Thyra::createMember(model->get_x_space());
305 Thyra::copy(*(model->getExactSolution(time).get_x()),solutionExact.ptr());
306 solutions.push_back(solutionExact);
307 auto solutionDotExact = Thyra::createMember(model->get_x_space());
308 Thyra::copy(*(model->getExactSolution(time).get_x_dot()),
309 solutionDotExact.ptr());
310 solutionsDot.push_back(solutionDotExact);
315 double xDotSlope = 0.0;
317 double order = stepper->getOrder();
320 solutions, xErrorNorm, xSlope,
321 solutionsDot, xDotErrorNorm, xDotSlope);
337 std::vector<RCP<Thyra::VectorBase<double>>> solutions;
338 std::vector<RCP<Thyra::VectorBase<double>>> solutionsDot;
339 std::vector<double> StepSize;
340 std::vector<double> xErrorNorm;
341 std::vector<double> xDotErrorNorm;
342 const int nTimeStepSizes = 4;
345 for (
int n=0; n<nTimeStepSizes; n++) {
349 getParametersFromXmlFile(
"Tempus_Trapezoidal_VanDerPol.xml");
357 if (n == nTimeStepSizes-1) dt /= 10.0;
362 .
sublist(
"Time Step Control").
set(
"Initial Time Step", dt);
363 integrator = Tempus::createIntegratorBasic<double>(pl, model);
366 bool integratorStatus = integrator->advanceTime();
370 time = integrator->getTime();
371 double timeFinal =pl->sublist(
"Demo Integrator")
372 .sublist(
"Time Step Control").
get<
double>(
"Final Time");
373 double tol = 100.0 * std::numeric_limits<double>::epsilon();
377 StepSize.push_back(dt);
378 auto solution = Thyra::createMember(model->get_x_space());
379 Thyra::copy(*(integrator->getX()),solution.ptr());
380 solutions.push_back(solution);
381 auto solutionDot = Thyra::createMember(model->get_x_space());
382 Thyra::copy(*(integrator->getXDot()),solutionDot.ptr());
383 solutionsDot.push_back(solutionDot);
387 if ((n == 0) || (n == nTimeStepSizes-1)) {
388 std::string fname =
"Tempus_Trapezoidal_VanDerPol-Ref.dat";
389 if (n == 0) fname =
"Tempus_Trapezoidal_VanDerPol.dat";
391 integrator->getSolutionHistory();
397 double xDotSlope = 0.0;
399 double order = stepper->getOrder();
402 solutions, xErrorNorm, xSlope,
403 solutionsDot, xDotErrorNorm, xDotSlope);
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)
#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)
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_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.