Tempus  Version of the Day
Time Integration
Tempus_ConvergenceTestUtils.hpp
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 #ifndef TEMPUS_CONVERGENCE_TEST_UTILS_HPP
10 #define TEMPUS_CONVERGENCE_TEST_UTILS_HPP
11 
12 #include <vector>
13 #include "Teuchos_as.hpp"
14 #include "Tempus_Stepper.hpp"
15 #include <fstream>
16 
17 
18 namespace Tempus_Test {
19 
20 /** \brief Linear regression class.
21  * Copied and modified from Rythmos.
22  */
23 template<class Scalar>
24 class LinearRegression
25 {
26  public:
27  LinearRegression();
28  void setData(std::vector<Scalar>& x, std::vector<Scalar>& y);
29  Scalar getSlope() const;
30  Scalar getYIntercept() const;
31  private:
32  // Private functions
33  void compute_();
34  void validateXYData_(std::vector<Scalar>& x, std::vector<Scalar>& y);
35 
36  // Private data
37  std::vector<Scalar> x_;
38  std::vector<Scalar> y_;
39  Scalar slope_;
40  Scalar yIntercept_;
41  bool isInitialized_;
42 };
43 
44 // Member functions:
45 
46 template<class Scalar>
47 LinearRegression<Scalar>::LinearRegression()
48 {
49  isInitialized_ = false;
50 }
51 
52 template<class Scalar>
53 void LinearRegression<Scalar>::setData(
54  std::vector<Scalar>& x, std::vector<Scalar>& y)
55 {
56  validateXYData_(x,y);
57  x_ = x; // copy x data
58  y_ = y; // copy y data
59  isInitialized_ = true;
60  compute_();
61 }
62 
63 template<class Scalar>
64 void LinearRegression<Scalar>::validateXYData_(
65  std::vector<Scalar>& x, std::vector<Scalar>& y)
66 {
67  TEUCHOS_TEST_FOR_EXCEPTION(x.size() != y.size(), std::logic_error,
68  "x and y data are note the same size for linear regression\n");
69  TEUCHOS_TEST_FOR_EXCEPTION(x.size() < 2, std::logic_error,
70  "Not enough data points for linear regression!\n");
71  int N = Teuchos::as<int>(x.size());
72  // There must be at least two unique x values
73  Scalar alpha = x[0];
74  int numUnique = 1;
75  for (int i=1; i<N ; ++i) {
76  if (x[i] != alpha) {
77  numUnique++;
78  }
79  }
80  TEUCHOS_TEST_FOR_EXCEPT(numUnique==1);
81 }
82 
83 template<class Scalar>
84 Scalar LinearRegression<Scalar>::getSlope() const
85 {
86  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
87  return slope_;
88 }
89 
90 template<class Scalar>
91 Scalar LinearRegression<Scalar>::getYIntercept() const
92 {
93  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
94  return yIntercept_;
95 }
96 
97 template<class Scalar>
98 void LinearRegression<Scalar>::compute_()
99 {
100  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
101  typedef Teuchos::ScalarTraits<Scalar> ST;
102 
103  int N = Teuchos::as<int>(x_.size());
104 
105  Scalar sum1 = ST::zero();
106  Scalar sum2 = ST::zero();
107  for (int i=0 ; i<N ; ++i) {
108  sum1 += x_[i]*y_[i];
109  sum2 += x_[i]*x_[i];
110  }
111  sum1 *= Scalar(-2*ST::one());
112  sum2 *= Scalar(-2*ST::one());
113 
114  Scalar sum3 = ST::zero();
115  Scalar sum4 = ST::zero();
116  for (int i=0 ; i<N ; ++i) {
117  for (int j=0 ; j<N ; ++j) {
118  sum3 += x_[i]*y_[j];
119  sum4 += x_[i]*x_[j];
120  }
121  }
122  sum3 *= Scalar(2*ST::one()/Scalar(N));
123  sum4 *= Scalar(2*ST::one()/Scalar(N));
124 
125  slope_ = ( sum3 + sum1 ) / ( sum4 + sum2 );
126 
127  yIntercept_ = ST::zero();
128  for (int i=0 ; i<N ; ++i ) {
129  yIntercept_ += y_[i]-slope_*x_[i];
130  }
131  yIntercept_ *= Scalar(ST::one()/Scalar(N));
132 }
133 
134 // Nonmember helper functions:
135 template<class Scalar>
136 Scalar computeLinearRegression(
137  std::vector<Scalar>& x, std::vector<Scalar>& y)
138 {
139  LinearRegression<Scalar> lr;
140  lr.setData(x,y);
141  return(lr.getSlope());
142 }
143 
144 template<class Scalar>
145 void computeLinearRegression(
146  std::vector<Scalar>& x, std::vector<Scalar>& y,
147  Scalar & slope, Scalar & yIntercept)
148 {
149  LinearRegression<Scalar> lr;
150  lr.setData(x,y);
151  slope = lr.getSlope();
152  yIntercept = lr.getYIntercept();
153  return;
154 }
155 
156 template<class Scalar>
157 Scalar computeLinearRegressionLogLog(
158  std::vector<Scalar>& x, std::vector<Scalar>& y)
159 {
160  TEUCHOS_TEST_FOR_EXCEPT(x.size() != y.size());
161  int N = Teuchos::as<int>(x.size());
162  std::vector<Scalar> xlog;
163  std::vector<Scalar> ylog;
164 
165  for (int i=0 ; i<N ; ++i) {
166  xlog.push_back(log(x[i]));
167  ylog.push_back(log(y[i]));
168  }
169 
170  LinearRegression<Scalar> lr;
171  lr.setData(xlog,ylog);
172  return(lr.getSlope());
173 }
174 
175 template<class Scalar>
176 Teuchos::RCP<LinearRegression<Scalar> > linearRegression()
177 {
178  Teuchos::RCP<LinearRegression<Scalar> > lr =
179  Teuchos::rcp(new LinearRegression<Scalar>());
180  return lr;
181 }
182 
183 template<class Scalar>
184 void writeOrderError(
185  const std::string filename,
186  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
187  std::vector<Scalar> & StepSize,
188  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
189  std::vector<Scalar> & xErrorNorm,
190  Scalar & xSlope,
191  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDot,
192  std::vector<Scalar> & xDotErrorNorm,
193  Scalar & xDotSlope,
194  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDotDot,
195  std::vector<Scalar> & xDotDotErrorNorm,
196  Scalar & xDotDotSlope)
197 {
198  Scalar order = stepper->getOrder();
199  int lastElem = solutions.size()-1; // Last element is the reference soln.
200  std::vector<Scalar> dt;
201 
202  auto ref_solution = solutions[lastElem];
203  for (std::size_t i=0; i < lastElem; ++i) {
204  dt.push_back(StepSize[i]);
205  auto tmp = solutions[i];
206  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
207  Scalar L2norm = Thyra::norm_2(*tmp);
208  xErrorNorm.push_back(L2norm);
209  }
210  xSlope = computeLinearRegressionLogLog<Scalar>(dt, xErrorNorm);
211 
212  if (!solutionsDot.empty()) {
213  auto ref_solutionDot = solutionsDot[lastElem];
214  for (std::size_t i=0; i < lastElem; ++i) {
215  auto tmp = solutionsDot[i];
216  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDot);
217  Scalar L2norm = Thyra::norm_2(*tmp);
218  xDotErrorNorm.push_back(L2norm);
219  }
220  xDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotErrorNorm);
221  }
222 
223  if (!solutionsDotDot.empty()) {
224  auto ref_solutionDotDot = solutionsDotDot[solutions.size()-1];
225  for (std::size_t i=0; i < lastElem; ++i) {
226  auto tmp = solutionsDotDot[i];
227  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDotDot);
228  Scalar L2norm = Thyra::norm_2(*tmp);
229  xDotDotErrorNorm.push_back(L2norm);
230  }
231  xDotDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotDotErrorNorm);
232  }
233 
234  std::cout << " Stepper = " << stepper->description() << std::endl;
235  std::cout << " =================================" << std::endl;
236  std::cout << " Expected Order = " << order << std::endl;
237  std::cout << " x Order = " << xSlope << std::endl;
238  if (!solutionsDot.empty())
239  std::cout<<" xDot Order = " << xDotSlope << std::endl;
240  if (!solutionsDotDot.empty())
241  std::cout<<" xDotDot Order = " << xDotDotSlope << std::endl;
242  std::cout << " =================================" << std::endl;
243 
244  std::ofstream ftmp(filename);
245  Scalar factor = 0.0;
246  for (int n=0; n<lastElem; n++) {
247  factor = 0.8*(pow(dt[n]/dt[0], order));
248  ftmp << dt[n]
249  << " " << xErrorNorm[n] << " " << factor*xErrorNorm[0];
250  if (xDotErrorNorm.size() == lastElem)
251  ftmp <<" "<< xDotErrorNorm[n] << " " << factor*xDotErrorNorm[0];
252  if (xDotDotErrorNorm.size() == lastElem)
253  ftmp <<" "<< xDotDotErrorNorm[n] << " " << factor*xDotDotErrorNorm[0];
254  ftmp << std::endl;
255  }
256  ftmp.close();
257 }
258 
259 template<class Scalar>
260 void writeOrderError(
261  const std::string filename,
262  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
263  std::vector<Scalar> & StepSize,
264  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
265  std::vector<Scalar> & xErrorNorm,
266  Scalar & xSlope,
267  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDot,
268  std::vector<Scalar> & xDotErrorNorm,
269  Scalar & xDotSlope)
270 {
271  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
272  std::vector<Scalar> xDotDotErrorNorm;
273  Scalar xDotDotSlope = Scalar(0.0);
274  writeOrderError(filename, stepper, StepSize,
275  solutions, xErrorNorm, xSlope,
276  solutionsDot, xDotErrorNorm, xDotSlope,
277  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
278 }
279 
280 template<class Scalar>
281 void writeOrderError(
282  const std::string filename,
283  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
284  std::vector<Scalar> & StepSize,
285  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
286  std::vector<Scalar> & xErrorNorm,
287  Scalar & xSlope)
288 {
289  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDot;
290  std::vector<Scalar> xDotErrorNorm;
291  Scalar xDotSlope = Scalar(0.0);
292  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
293  std::vector<Scalar> xDotDotErrorNorm;
294  Scalar xDotDotSlope = Scalar(0.0);
295  writeOrderError(filename, stepper, StepSize,
296  solutions, xErrorNorm, xSlope,
297  solutionsDot, xDotErrorNorm, xDotSlope,
298  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
299 }
300 
301 template<class Scalar>
302 void writeSolution(
303  const std::string filename,
304  Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > solutionHistory)
305 {
306  std::ofstream ftmp(filename);
307  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x;
308  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDot;
309  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDotDot;
310  for (int i=0; i<solutionHistory->getNumStates(); i++) {
311  Teuchos::RCP<const Tempus::SolutionState<Scalar> > solutionState =
312  (*solutionHistory)[i];
313  Scalar time = solutionState->getTime();
314  x = solutionState->getX();
315  xDot = solutionState->getXDot();
316  xDotDot = solutionState->getXDotDot();
317  int J = x->space()->dim();
318 
319  ftmp << time;
320  for (int j=0; j<J; j++) ftmp << " " << get_ele(*x, j);
321  if (xDot != Teuchos::null)
322  for (int j=0; j<J; j++) ftmp << " " << get_ele(*xDot, j);
323  if (xDotDot != Teuchos::null)
324  for (int j=0; j<J; j++) ftmp << " " << get_ele(*xDotDot, j);
325  ftmp << std::endl;
326  }
327  ftmp.close();
328 }
329 
330 template<class Scalar>
331 void writeSolution(
332  const std::string filename,
333  Teuchos::RCP<Tempus::SolutionHistory<Scalar> > solutionHistory)
334 {
335  Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > sh(solutionHistory);
336  writeSolution(filename, sh);
337 }
338 
339 } // namespace Tempus_Test
340 
341 #endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP
342 
Tempus_Test::writeSolution
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Definition: Tempus_ConvergenceTestUtils.hpp:302
Tempus_Test::writeOrderError
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)
Definition: Tempus_ConvergenceTestUtils.hpp:184
Tempus_Test
Definition: Tempus_BackwardEuler_ASA.cpp:34
Tempus_Test::computeLinearRegression
Scalar computeLinearRegression(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition: Tempus_ConvergenceTestUtils.hpp:136
Tempus_Test::LinearRegression::compute_
void compute_()
Definition: Tempus_ConvergenceTestUtils.hpp:98
Tempus::solutionHistory
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
Definition: Tempus_SolutionHistory_impl.hpp:504
Tempus_Test::LinearRegression::validateXYData_
void validateXYData_(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition: Tempus_ConvergenceTestUtils.hpp:64
Tempus_Test::LinearRegression::x_
std::vector< Scalar > x_
Definition: Tempus_ConvergenceTestUtils.hpp:37
Tempus_Test::LinearRegression::yIntercept_
Scalar yIntercept_
Definition: Tempus_ConvergenceTestUtils.hpp:40
Tempus_Test::LinearRegression
Linear regression class. Copied and modified from Rythmos.
Definition: Tempus_ConvergenceTestUtils.hpp:24
Tempus_Test::LinearRegression::LinearRegression
LinearRegression()
Definition: Tempus_ConvergenceTestUtils.hpp:47
Tempus_Test::LinearRegression::slope_
Scalar slope_
Definition: Tempus_ConvergenceTestUtils.hpp:39
Tempus_Stepper.hpp
Tempus_Test::LinearRegression::getYIntercept
Scalar getYIntercept() const
Definition: Tempus_ConvergenceTestUtils.hpp:91
Tempus_Test::LinearRegression::setData
void setData(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition: Tempus_ConvergenceTestUtils.hpp:53
Tempus_Test::LinearRegression::isInitialized_
bool isInitialized_
Definition: Tempus_ConvergenceTestUtils.hpp:41
Tempus_Test::LinearRegression::getSlope
Scalar getSlope() const
Definition: Tempus_ConvergenceTestUtils.hpp:84
Tempus::SolutionHistory
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Definition: Tempus_Integrator.hpp:25
Tempus::Stepper
Thyra Base interface for time steppers.
Definition: Tempus_Integrator.hpp:24
Tempus_Test::LinearRegression::y_
std::vector< Scalar > y_
Definition: Tempus_ConvergenceTestUtils.hpp:38
Tempus_Test::linearRegression
Teuchos::RCP< LinearRegression< Scalar > > linearRegression()
Definition: Tempus_ConvergenceTestUtils.hpp:176
Tempus_Test::computeLinearRegressionLogLog
Scalar computeLinearRegressionLogLog(std::vector< Scalar > &x, std::vector< Scalar > &y)
Definition: Tempus_ConvergenceTestUtils.hpp:157