Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_PhysicsStateTest.cpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#include "Teuchos_UnitTestHarness.hpp"
10#include "Teuchos_XMLParameterListHelpers.hpp"
11#include "Teuchos_TimeMonitor.hpp"
12
13#include "Thyra_VectorStdOps.hpp"
14
15#include "Tempus_IntegratorBasic.hpp"
18
19#include "../TestModels/SinCosModel.hpp"
20#include "../TestModels/VanDerPolModel.hpp"
22
23#include <fstream>
24#include <vector>
25
26namespace Tempus_Test {
27
28using Teuchos::RCP;
29using Teuchos::ParameterList;
30using Teuchos::sublist;
31using Teuchos::getParametersFromXmlFile;
32
36
37
38// ************************************************************
39// ************************************************************
40TEUCHOS_UNIT_TEST(PhysicsState, SinCos)
41{
42 std::vector<double> StepSize;
43 std::vector<double> ErrorNorm;
44 const int nTimeStepSizes = 7;
45 double dt = 0.2;
46 double order = 0.0;
47 for (int n=0; n<nTimeStepSizes; n++) {
48
49 // Read params from .xml file
50 RCP<ParameterList> pList =
51 getParametersFromXmlFile("Tempus_PhysicsState_SinCos.xml");
52
53 //std::ofstream ftmp("PL.txt");
54 //pList->print(ftmp);
55 //ftmp.close();
56
57 // Setup the SinCosModel
58 RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
59 //RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
60 RCP<SinCosModel<double> > model =
61 Teuchos::rcp(new SinCosModel<double> (scm_pl));
62
63 dt /= 2;
64
65 // Setup the Integrator and reset initial time step
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);
71
72 // Replace Tempus::StepperForwardEuler with
73 // Tempus_Test::StepperPhysicsStateTest
74 Teuchos::RCP<Tempus::Stepper<double> > physicsStepper = Teuchos::rcp(
75 new StepperPhysicsStateTest<double>(model));
76 integrator->setStepper(physicsStepper);
77 order = integrator->getStepper()->getOrder();
78
79 // Initial Conditions
80 // During the Integrator construction, the initial SolutionState
81 // is set by default to model->getNominalVales().get_x(). However,
82 // the application can set it also by integrator->initializeSolutionHistory.
83 RCP<Thyra::VectorBase<double> > x0 =
84 model->getNominalValues().get_x()->clone_v();
85 integrator->initializeSolutionHistory(0.0, x0);
86
87
88 // Replace Tempus::PhysicsState with
89 // Tempus_Test::PhysicsStateCounter
90 RCP<PhysicsStateCounter<double> > pSC = Teuchos::rcp(
91 new PhysicsStateCounter<double> ("PhysicsStateTest", 0));
92 integrator->getSolutionHistory()->getCurrentState()->setPhysicsState(pSC);
93
94 // Integrate to timeMax
95 bool integratorStatus = integrator->advanceTime();
96 TEST_ASSERT(integratorStatus)
97
98 // Test PhysicsState
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);
103 //out << " Counter = " << pSC->getCounter() << std::endl;
104 TEST_EQUALITY(pSC->getCounter(), (int)(1.0/dt));
105
106
107 // Test if at 'Final Time'
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);
112
113 // Time-integrated solution and the exact solution
114 RCP<Thyra::VectorBase<double> > x = integrator->getX();
115 RCP<const Thyra::VectorBase<double> > x_exact =
116 model->getExactSolution(time).get_x();
117
118 // Plot sample solution and exact solution
119 if (n == 0) {
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;
134 }
135 ftmp.close();
136 }
137
138 // Calculate the error
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);
144 }
145
146 // Check the order and intercept
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 );
155
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;
161 }
162 ftmp.close();
163
164 Teuchos::TimeMonitor::summarize();
165}
166
167} // namespace Tempus_Test
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
TEUCHOS_UNIT_TEST(BackwardEuler, SinCos_ASA)