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 namespace Tempus_Test {
21 
25 template <class Scalar>
27  public:
29  void setData(std::vector<Scalar>& x, std::vector<Scalar>& y);
30  Scalar getSlope() const;
31  Scalar getYIntercept() const;
32 
33  private:
34  // Private functions
35  void compute_();
36  void validateXYData_(std::vector<Scalar>& x, std::vector<Scalar>& y);
37 
38  // Private data
39  std::vector<Scalar> x_;
40  std::vector<Scalar> y_;
41  Scalar slope_;
42  Scalar yIntercept_;
44 };
45 
46 // Member functions:
47 
48 template <class Scalar>
50 {
51  isInitialized_ = false;
52 }
53 
54 template <class Scalar>
55 void LinearRegression<Scalar>::setData(std::vector<Scalar>& x,
56  std::vector<Scalar>& y)
57 {
58  validateXYData_(x, y);
59  x_ = x; // copy x data
60  y_ = y; // copy y data
61  isInitialized_ = true;
62  compute_();
63 }
64 
65 template <class Scalar>
66 void LinearRegression<Scalar>::validateXYData_(std::vector<Scalar>& x,
67  std::vector<Scalar>& y)
68 {
70  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>
139 Scalar computeLinearRegression(std::vector<Scalar>& x, std::vector<Scalar>& y)
140 {
142  lr.setData(x, y);
143  return (lr.getSlope());
144 }
145 
146 template <class Scalar>
147 void computeLinearRegression(std::vector<Scalar>& x, std::vector<Scalar>& y,
148  Scalar& slope, Scalar& yIntercept)
149 {
151  lr.setData(x, y);
152  slope = lr.getSlope();
153  yIntercept = lr.getYIntercept();
154  return;
155 }
156 
157 template <class Scalar>
158 Scalar computeLinearRegressionLogLog(std::vector<Scalar>& x,
159  std::vector<Scalar>& y)
160 {
161  TEUCHOS_TEST_FOR_EXCEPT(x.size() != y.size());
162  int N = Teuchos::as<int>(x.size());
163  std::vector<Scalar> xlog;
164  std::vector<Scalar> ylog;
165 
166  for (int i = 0; i < N; ++i) {
167  if (!(Tempus::approxZero(x[i]) || Tempus::approxZero(y[i]))) {
168  xlog.push_back(log(x[i]));
169  ylog.push_back(log(y[i]));
170  }
171  }
172 
174  lr.setData(xlog, ylog);
175  return (lr.getSlope());
176 }
177 
178 template <class Scalar>
180 {
183  return lr;
184 }
185 
186 template <class Scalar>
188  const std::string filename, Teuchos::RCP<Tempus::Stepper<Scalar>> stepper,
189  std::vector<Scalar>& StepSize,
190  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>>& solutions,
191  std::vector<Scalar>& xErrorNorm, Scalar& xSlope,
192  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>>& solutionsDot,
193  std::vector<Scalar>& xDotErrorNorm, Scalar& xDotSlope,
194  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>>& solutionsDotDot,
195  std::vector<Scalar>& xDotDotErrorNorm, Scalar& xDotDotSlope,
197 {
198  if (solutions.empty()) {
199  out << "No solutions to write out!\n"
200  << std::endl;
201  return;
202  }
203  std::size_t lastElem = solutions.size() - 1; // Last element is the ref soln.
204  std::vector<Scalar> dt;
205 
206  auto ref_solution = solutions[lastElem];
207  for (std::size_t i = 0; i < lastElem; ++i) {
208  dt.push_back(StepSize[i]);
209  auto tmp = solutions[i];
210  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solution);
211  Scalar L2norm = Thyra::norm_2(*tmp);
212  xErrorNorm.push_back(L2norm);
213  }
214  xSlope = computeLinearRegressionLogLog<Scalar>(dt, xErrorNorm);
215 
216  if (!solutionsDot.empty()) {
217  auto ref_solutionDot = solutionsDot[lastElem];
218  for (std::size_t i = 0; i < lastElem; ++i) {
219  auto tmp = solutionsDot[i];
220  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDot);
221  Scalar L2norm = Thyra::norm_2(*tmp);
222  xDotErrorNorm.push_back(L2norm);
223  }
224  xDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotErrorNorm);
225  }
226 
227  if (!solutionsDotDot.empty()) {
228  auto ref_solutionDotDot = solutionsDotDot[solutions.size() - 1];
229  for (std::size_t i = 0; i < lastElem; ++i) {
230  auto tmp = solutionsDotDot[i];
231  Thyra::Vp_StV(tmp.ptr(), -1.0, *ref_solutionDotDot);
232  Scalar L2norm = Thyra::norm_2(*tmp);
233  xDotDotErrorNorm.push_back(L2norm);
234  }
235  xDotDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotDotErrorNorm);
236  }
237 
238  Scalar order = stepper->getOrder();
239  out << " Stepper = " << stepper->description() << std::endl;
240  out << " =================================" << std::endl;
241  out << " Expected Order = " << order << std::endl;
242  out << " x Order = " << xSlope << std::endl;
243  if (!solutionsDot.empty())
244  out << " xDot Order = " << xDotSlope << std::endl;
245  if (!solutionsDotDot.empty())
246  out << " xDotDot Order = " << xDotDotSlope << std::endl;
247  out << " =================================" << std::endl;
248 
249  std::ofstream ftmp(filename);
250  Scalar factor = 0.0;
251  for (std::size_t n = 0; n < lastElem; n++) {
252  factor = 0.8 * (pow(dt[n] / dt[0], order));
253  ftmp << dt[n] << " " << 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] << " "
258  << factor * xDotDotErrorNorm[0];
259  ftmp << std::endl;
260  }
261  ftmp.close();
262 }
263 
264 template <class Scalar>
266  const std::string filename, 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, Scalar& xSlope,
270  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>>& solutionsDot,
271  std::vector<Scalar>& xDotErrorNorm, Scalar& xDotSlope,
273 {
274  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>> solutionsDotDot;
275  std::vector<Scalar> xDotDotErrorNorm;
276  Scalar xDotDotSlope = Scalar(0.0);
277  writeOrderError(filename, stepper, StepSize, solutions, xErrorNorm, xSlope,
278  solutionsDot, xDotErrorNorm, xDotSlope, solutionsDotDot,
279  xDotDotErrorNorm, xDotDotSlope, out);
280 }
281 
282 template <class Scalar>
284  const std::string filename, Teuchos::RCP<Tempus::Stepper<Scalar>> stepper,
285  std::vector<Scalar>& StepSize,
286  std::vector<Teuchos::RCP<Thyra::VectorBase<Scalar>>>& solutions,
287  std::vector<Scalar>& xErrorNorm, Scalar& xSlope, Teuchos::FancyOStream& out)
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, solutions, xErrorNorm, xSlope,
296  solutionsDot, xDotErrorNorm, xDotSlope, solutionsDotDot,
297  xDotDotErrorNorm, xDotDotSlope, out);
298 }
299 
300 template <class Scalar>
302  const std::string filename,
303  Teuchos::RCP<const Tempus::SolutionHistory<Scalar>> solutionHistory)
304 {
305  std::ofstream ftmp(filename);
309  for (int i = 0; i < solutionHistory->getNumStates(); i++) {
311  (*solutionHistory)[i];
312  Scalar time = solutionState->getTime();
313  x = solutionState->getX();
314  xDot = solutionState->getXDot();
315  xDotDot = solutionState->getXDotDot();
316  int J = x->space()->dim();
317 
318  ftmp << time;
319  for (int j = 0; j < J; j++) ftmp << " " << get_ele(*x, j);
320  if (xDot != Teuchos::null)
321  for (int j = 0; j < J; j++) ftmp << " " << get_ele(*xDot, j);
322  if (xDotDot != Teuchos::null)
323  for (int j = 0; j < J; j++) ftmp << " " << get_ele(*xDotDot, j);
324  ftmp << std::endl;
325  }
326  ftmp.close();
327 }
328 
329 template <class Scalar>
330 void writeSolution(
331  const std::string filename,
333 {
335  writeSolution(filename, sh);
336 }
337 
338 } // namespace Tempus_Test
339 
340 #endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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::FancyOStream &out)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar >> solutionHistory)
Thyra Base interface for time steppers.
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)