42 std::vector<double> StepSize;
43 std::vector<double> ErrorNorm;
44 const int nTimeStepSizes = 7;
47 for (
int n=0; n<nTimeStepSizes; n++) {
50 RCP<ParameterList> pList =
51 getParametersFromXmlFile(
"Tempus_PhysicsState_SinCos.xml");
58 RCP<ParameterList> scm_pl = sublist(pList,
"SinCosModel",
true);
60 RCP<SinCosModel<double> > model =
61 Teuchos::rcp(
new SinCosModel<double> (scm_pl));
66 RCP<ParameterList> pl = sublist(pList,
"Tempus",
true);
67 pl->sublist(
"Demo Integrator")
68 .sublist(
"Time Step Control").set(
"Initial Time Step", dt);
69 RCP<Tempus::IntegratorBasic<double> > integrator =
70 Tempus::createIntegratorBasic<double>(pl, model);
74 Teuchos::RCP<Tempus::Stepper<double> > physicsStepper = Teuchos::rcp(
75 new StepperPhysicsStateTest<double>(model));
76 integrator->setStepper(physicsStepper);
77 order = integrator->getStepper()->getOrder();
83 RCP<Thyra::VectorBase<double> > x0 =
84 model->getNominalValues().get_x()->clone_v();
85 integrator->initializeSolutionHistory(0.0, x0);
90 RCP<PhysicsStateCounter<double> > pSC = Teuchos::rcp(
91 new PhysicsStateCounter<double> (
"PhysicsStateTest", 0));
92 integrator->getSolutionHistory()->getCurrentState()->setPhysicsState(pSC);
95 bool integratorStatus = integrator->advanceTime();
96 TEST_ASSERT(integratorStatus)
99 Teuchos::RCP<Tempus::PhysicsState<double> > pS =
100 integrator->getSolutionHistory()->getCurrentState()->getPhysicsState();
101 TEST_EQUALITY(pS->getName(),
"PhysicsStateTest");
102 pSC = Teuchos::rcp_dynamic_cast<PhysicsStateCounter<double> >(pS);
104 TEST_EQUALITY(pSC->getCounter(), (
int)(1.0/dt));
108 double time = integrator->getTime();
109 double timeFinal = pl->sublist(
"Demo Integrator")
110 .sublist(
"Time Step Control").get<
double>(
"Final Time");
111 TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
114 RCP<Thyra::VectorBase<double> > x = integrator->getX();
115 RCP<const Thyra::VectorBase<double> > x_exact =
116 model->getExactSolution(time).get_x();
120 std::ofstream ftmp(
"Tempus_ForwardEuler_SinCos.dat");
121 RCP<const SolutionHistory<double> > solutionHistory =
122 integrator->getSolutionHistory();
123 RCP<const Thyra::VectorBase<double> > x_exact_plot;
124 for (
int i=0; i<solutionHistory->getNumStates(); i++) {
125 RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
126 double time_i = solutionState->getTime();
127 RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
128 x_exact_plot = model->getExactSolution(time_i).get_x();
129 ftmp << time_i <<
" "
130 << Thyra::get_ele(*(x_plot), 0) <<
" "
131 << Thyra::get_ele(*(x_plot), 1) <<
" "
132 << Thyra::get_ele(*(x_exact_plot), 0) <<
" "
133 << Thyra::get_ele(*(x_exact_plot), 1) << std::endl;
139 RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
140 Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
141 StepSize.push_back(dt);
142 const double L2norm = Thyra::norm_2(*xdiff);
143 ErrorNorm.push_back(L2norm);
147 double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
148 out <<
" Stepper = ForwardEuler" << std::endl;
149 out <<
" =========================" << std::endl;
150 out <<
" Expected order: " << order << std::endl;
151 out <<
" Observed order: " << slope << std::endl;
152 out <<
" =========================" << std::endl;
153 TEST_FLOATING_EQUALITY( slope, order, 0.01 );
154 TEST_FLOATING_EQUALITY( ErrorNorm[0], 0.051123, 1.0e-4 );
156 std::ofstream ftmp(
"Tempus_ForwardEuler_SinCos-Error.dat");
157 double error0 = 0.8*ErrorNorm[0];
158 for (
int n=0; n<nTimeStepSizes; n++) {
159 ftmp << StepSize[n] <<
" " << ErrorNorm[n] <<
" "
160 << error0*(pow(StepSize[n]/StepSize[0],order)) << std::endl;
164 Teuchos::TimeMonitor::summarize();
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...