Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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"
21 #include "../TestUtils/Tempus_ConvergenceTestUtils.hpp"
22 
23 #include <fstream>
24 #include <vector>
25 
26 namespace Tempus_Test {
27 
28 using Teuchos::getParametersFromXmlFile;
30 using Teuchos::RCP;
31 using Teuchos::sublist;
32 
36 
37 // ************************************************************
38 // ************************************************************
39 TEUCHOS_UNIT_TEST(PhysicsState, SinCos)
40 {
41  std::vector<double> StepSize;
42  std::vector<double> ErrorNorm;
43  const int nTimeStepSizes = 7;
44  double dt = 0.2;
45  double order = 0.0;
46  for (int n = 0; n < nTimeStepSizes; n++) {
47  // Read params from .xml file
48  RCP<ParameterList> pList =
49  getParametersFromXmlFile("Tempus_PhysicsState_SinCos.xml");
50 
51  // std::ofstream ftmp("PL.txt");
52  // pList->print(ftmp);
53  // ftmp.close();
54 
55  // Setup the SinCosModel
56  RCP<ParameterList> scm_pl = sublist(pList, "SinCosModel", true);
57  // RCP<SinCosModel<double> > model = sineCosineModel(scm_pl);
58  RCP<SinCosModel<double> > model =
59  Teuchos::rcp(new SinCosModel<double>(scm_pl));
60 
61  dt /= 2;
62 
63  // Setup the Integrator and reset initial time step
64  RCP<ParameterList> pl = sublist(pList, "Tempus", true);
65  pl->sublist("Demo Integrator")
66  .sublist("Time Step Control")
67  .set("Initial Time Step", dt);
69  Tempus::createIntegratorBasic<double>(pl, model);
70 
71  // Replace Tempus::StepperForwardEuler with
72  // Tempus_Test::StepperPhysicsStateTest
73  Teuchos::RCP<Tempus::Stepper<double> > physicsStepper =
75  integrator->setStepper(physicsStepper);
76  order = integrator->getStepper()->getOrder();
77 
78  // Initial Conditions
79  // During the Integrator construction, the initial SolutionState
80  // is set by default to model->getNominalVales().get_x(). However,
81  // the application can set it also by integrator->initializeSolutionHistory.
83  model->getNominalValues().get_x()->clone_v();
84  integrator->initializeSolutionHistory(0.0, x0);
85 
86  // Replace Tempus::PhysicsState with
87  // Tempus_Test::PhysicsStateCounter
89  Teuchos::rcp(new PhysicsStateCounter<double>("PhysicsStateTest", 0));
90  integrator->getSolutionHistory()->getCurrentState()->setPhysicsState(pSC);
91 
92  // Integrate to timeMax
93  bool integratorStatus = integrator->advanceTime();
94  TEST_ASSERT(integratorStatus)
95 
96  // Test PhysicsState
98  integrator->getSolutionHistory()->getCurrentState()->getPhysicsState();
99  TEST_EQUALITY(pS->getName(), "PhysicsStateTest");
100  pSC = Teuchos::rcp_dynamic_cast<PhysicsStateCounter<double> >(pS);
101  // out << " Counter = " << pSC->getCounter() << std::endl;
102  TEST_EQUALITY(pSC->getCounter(), (int)(1.0 / dt));
103 
104  // Test if at 'Final Time'
105  double time = integrator->getTime();
106  double timeFinal = pl->sublist("Demo Integrator")
107  .sublist("Time Step Control")
108  .get<double>("Final Time");
109  TEST_FLOATING_EQUALITY(time, timeFinal, 1.0e-14);
110 
111  // Time-integrated solution and the exact solution
112  RCP<Thyra::VectorBase<double> > x = integrator->getX();
114  model->getExactSolution(time).get_x();
115 
116  // Plot sample solution and exact solution
117  if (n == 0) {
118  std::ofstream ftmp("Tempus_ForwardEuler_SinCos.dat");
119  RCP<const SolutionHistory<double> > solutionHistory =
120  integrator->getSolutionHistory();
122  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
123  RCP<const SolutionState<double> > solutionState = (*solutionHistory)[i];
124  double time_i = solutionState->getTime();
125  RCP<const Thyra::VectorBase<double> > x_plot = solutionState->getX();
126  x_exact_plot = model->getExactSolution(time_i).get_x();
127  ftmp << time_i << " " << Thyra::get_ele(*(x_plot), 0) << " "
128  << Thyra::get_ele(*(x_plot), 1) << " "
129  << Thyra::get_ele(*(x_exact_plot), 0) << " "
130  << Thyra::get_ele(*(x_exact_plot), 1) << std::endl;
131  }
132  ftmp.close();
133  }
134 
135  // Calculate the error
136  RCP<Thyra::VectorBase<double> > xdiff = x->clone_v();
137  Thyra::V_StVpStV(xdiff.ptr(), 1.0, *x_exact, -1.0, *(x));
138  StepSize.push_back(dt);
139  const double L2norm = Thyra::norm_2(*xdiff);
140  ErrorNorm.push_back(L2norm);
141  }
142 
143  // Check the order and intercept
144  double slope = computeLinearRegressionLogLog<double>(StepSize, ErrorNorm);
145  out << " Stepper = ForwardEuler" << std::endl;
146  out << " =========================" << std::endl;
147  out << " Expected order: " << order << std::endl;
148  out << " Observed order: " << slope << std::endl;
149  out << " =========================" << std::endl;
150  TEST_FLOATING_EQUALITY(slope, order, 0.01);
151  TEST_FLOATING_EQUALITY(ErrorNorm[0], 0.051123, 1.0e-4);
152 
153  std::ofstream ftmp("Tempus_ForwardEuler_SinCos-Error.dat");
154  double error0 = 0.8 * ErrorNorm[0];
155  for (int n = 0; n < nTimeStepSizes; n++) {
156  ftmp << StepSize[n] << " " << ErrorNorm[n] << " "
157  << error0 * (pow(StepSize[n] / StepSize[0], order)) << std::endl;
158  }
159  ftmp.close();
160 
162 }
163 
164 } // namespace Tempus_Test
This is a Forward Euler time stepper to test the PhysicsState.
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)
#define TEST_ASSERT(v1)
T * get() const
Sine-Cosine model problem from Rythmos. This is a canonical Sine-Cosine differential equation with a...
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)
Ptr< T > ptr() const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
PhysicsStateCounter is a simple PhysicsState that counts steps.
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEST_EQUALITY(v1, v2)
Solution state for integrators and steppers.