Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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>
25 {
26  public:
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_;
42 };
43 
44 // Member functions:
45 
46 template<class Scalar>
48 {
49  isInitialized_ = false;
50 }
51 
52 template<class Scalar>
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>
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>
85 {
86  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
87  return slope_;
88 }
89 
90 template<class Scalar>
92 {
93  TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
94  return yIntercept_;
95 }
96 
97 template<class Scalar>
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>
137  std::vector<Scalar>& x, std::vector<Scalar>& y)
138 {
140  lr.setData(x,y);
141  return(lr.getSlope());
142 }
143 
144 template<class Scalar>
146  std::vector<Scalar>& x, std::vector<Scalar>& y,
147  Scalar & slope, Scalar & yIntercept)
148 {
150  lr.setData(x,y);
151  slope = lr.getSlope();
152  yIntercept = lr.getYIntercept();
153  return;
154 }
155 
156 template<class Scalar>
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 
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>
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  if ( solutions.empty() ) {
199  std::cout << "No solutions to write out!\n" << std::endl;
200  return;
201  }
202  std::size_t lastElem = solutions.size()-1; // Last element is the ref soln.
203  std::vector<Scalar> dt;
204 
205  auto ref_solution = solutions[lastElem];
206  for (std::size_t i=0; i < lastElem; ++i) {
207  dt.push_back(StepSize[i]);
208  auto tmp = solutions[i];
209  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
210  Scalar L2norm = Thyra::norm_2(*tmp);
211  xErrorNorm.push_back(L2norm);
212  }
213  xSlope = computeLinearRegressionLogLog<Scalar>(dt, xErrorNorm);
214 
215  if (!solutionsDot.empty()) {
216  auto ref_solutionDot = solutionsDot[lastElem];
217  for (std::size_t i=0; i < lastElem; ++i) {
218  auto tmp = solutionsDot[i];
219  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDot);
220  Scalar L2norm = Thyra::norm_2(*tmp);
221  xDotErrorNorm.push_back(L2norm);
222  }
223  xDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotErrorNorm);
224  }
225 
226  if (!solutionsDotDot.empty()) {
227  auto ref_solutionDotDot = solutionsDotDot[solutions.size()-1];
228  for (std::size_t i=0; i < lastElem; ++i) {
229  auto tmp = solutionsDotDot[i];
230  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDotDot);
231  Scalar L2norm = Thyra::norm_2(*tmp);
232  xDotDotErrorNorm.push_back(L2norm);
233  }
234  xDotDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotDotErrorNorm);
235  }
236 
237  Scalar order = stepper->getOrder();
238  std::cout << " Stepper = " << stepper->description() << std::endl;
239  std::cout << " =================================" << std::endl;
240  std::cout << " Expected Order = " << order << std::endl;
241  std::cout << " x Order = " << xSlope << std::endl;
242  if (!solutionsDot.empty())
243  std::cout<<" xDot Order = " << xDotSlope << std::endl;
244  if (!solutionsDotDot.empty())
245  std::cout<<" xDotDot Order = " << xDotDotSlope << std::endl;
246  std::cout << " =================================" << std::endl;
247 
248  std::ofstream ftmp(filename);
249  Scalar factor = 0.0;
250  for (std::size_t n=0; n<lastElem; n++) {
251  factor = 0.8*(pow(dt[n]/dt[0], order));
252  ftmp << dt[n]
253  << " " << xErrorNorm[n] << " " << factor*xErrorNorm[0];
254  if (xDotErrorNorm.size() == lastElem)
255  ftmp <<" "<< xDotErrorNorm[n] << " " << factor*xDotErrorNorm[0];
256  if (xDotDotErrorNorm.size() == lastElem)
257  ftmp <<" "<< xDotDotErrorNorm[n] << " " << factor*xDotDotErrorNorm[0];
258  ftmp << std::endl;
259  }
260  ftmp.close();
261 }
262 
263 template<class Scalar>
265  const std::string filename,
266  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
267  std::vector<Scalar> & StepSize,
268  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
269  std::vector<Scalar> & xErrorNorm,
270  Scalar & xSlope,
271  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutionsDot,
272  std::vector<Scalar> & xDotErrorNorm,
273  Scalar & xDotSlope)
274 {
275  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
276  std::vector<Scalar> xDotDotErrorNorm;
277  Scalar xDotDotSlope = Scalar(0.0);
278  writeOrderError(filename, stepper, StepSize,
279  solutions, xErrorNorm, xSlope,
280  solutionsDot, xDotErrorNorm, xDotSlope,
281  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
282 }
283 
284 template<class Scalar>
286  const std::string filename,
287  Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
288  std::vector<Scalar> & StepSize,
289  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> & solutions,
290  std::vector<Scalar> & xErrorNorm,
291  Scalar & xSlope)
292 {
293  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDot;
294  std::vector<Scalar> xDotErrorNorm;
295  Scalar xDotSlope = Scalar(0.0);
296  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
297  std::vector<Scalar> xDotDotErrorNorm;
298  Scalar xDotDotSlope = Scalar(0.0);
299  writeOrderError(filename, stepper, StepSize,
300  solutions, xErrorNorm, xSlope,
301  solutionsDot, xDotErrorNorm, xDotSlope,
302  solutionsDotDot, xDotDotErrorNorm, xDotDotSlope);
303 }
304 
305 template<class Scalar>
307  const std::string filename,
309 {
310  std::ofstream ftmp(filename);
311  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x;
312  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDot;
313  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDotDot;
314  for (int i=0; i<solutionHistory->getNumStates(); i++) {
315  Teuchos::RCP<const Tempus::SolutionState<Scalar> > solutionState =
316  (*solutionHistory)[i];
317  Scalar time = solutionState->getTime();
318  x = solutionState->getX();
319  xDot = solutionState->getXDot();
320  xDotDot = solutionState->getXDotDot();
321  int J = x->space()->dim();
322 
323  ftmp << time;
324  for (int j=0; j<J; j++) ftmp << " " << get_ele(*x, j);
325  if (xDot != Teuchos::null)
326  for (int j=0; j<J; j++) ftmp << " " << get_ele(*xDot, j);
327  if (xDotDot != Teuchos::null)
328  for (int j=0; j<J; j++) ftmp << " " << get_ele(*xDotDot, j);
329  ftmp << std::endl;
330  }
331  ftmp.close();
332 }
333 
334 template<class Scalar>
335 void writeSolution(
336  const std::string filename,
338 {
339  Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > sh(solutionHistory);
340  writeSolution(filename, sh);
341 }
342 
343 } // namespace Tempus_Test
344 
345 #endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP
346 
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Thyra Base interface for time steppers.
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)
Scalar computeLinearRegressionLogLog(std::vector< Scalar > &x, std::vector< Scalar > &y)
Teuchos::RCP< LinearRegression< Scalar > > linearRegression()
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
void validateXYData_(std::vector< Scalar > &x, std::vector< Scalar > &y)
Scalar computeLinearRegression(std::vector< Scalar > &x, std::vector< Scalar > &y)
void setData(std::vector< Scalar > &x, std::vector< Scalar > &y)
Linear regression class. Copied and modified from Rythmos.