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