Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
21namespace Tempus_Test {
22
26template<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_;
45};
46
47// Member functions:
48
49template<class Scalar>
51{
52 isInitialized_ = false;
53}
54
55template<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
66template<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
86template<class Scalar>
88{
89 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
90 return slope_;
91}
92
93template<class Scalar>
95{
96 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
97 return yIntercept_;
98}
99
100template<class Scalar>
102{
103 TEUCHOS_TEST_FOR_EXCEPT(!isInitialized_);
104 typedef Teuchos::ScalarTraits<Scalar> ST;
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:
138template<class Scalar>
140 std::vector<Scalar>& x, std::vector<Scalar>& y)
141{
143 lr.setData(x,y);
144 return(lr.getSlope());
145}
146
147template<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
159template<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
180template<class Scalar>
181Teuchos::RCP<LinearRegression<Scalar> > linearRegression()
182{
183 Teuchos::RCP<LinearRegression<Scalar> > lr =
184 Teuchos::rcp(new LinearRegression<Scalar>());
185 return lr;
186}
187
188template<class Scalar>
190 const std::string filename,
191 Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
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
268template<class Scalar>
270 const std::string filename,
271 Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
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
289template<class Scalar>
291 const std::string filename,
292 Teuchos::RCP<Tempus::Stepper<Scalar> > stepper,
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
310template<class Scalar>
312 const std::string filename,
313 Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > solutionHistory)
314{
315 std::ofstream ftmp(filename);
316 Teuchos::RCP<const Thyra::VectorBase<Scalar> > x;
317 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDot;
318 Teuchos::RCP<const Thyra::VectorBase<Scalar> > xDotDot;
319 for (int i=0; i<solutionHistory->getNumStates(); i++) {
320 Teuchos::RCP<const Tempus::SolutionState<Scalar> > solutionState =
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
339template<class Scalar>
340void writeSolution(
341 const std::string filename,
342 Teuchos::RCP<Tempus::SolutionHistory<Scalar> > solutionHistory)
343{
344 Teuchos::RCP<const Tempus::SolutionHistory<Scalar> > sh(solutionHistory);
345 writeSolution(filename, sh);
346}
347
348} // namespace Tempus_Test
349
350#endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP
351
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
Thyra Base interface for time steppers.
Linear regression class. Copied and modified from Rythmos.
void validateXYData_(std::vector< Scalar > &x, std::vector< Scalar > &y)
void setData(std::vector< Scalar > &x, std::vector< Scalar > &y)
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)
void writeSolution(const std::string filename, Teuchos::RCP< const Tempus::SolutionHistory< Scalar > > solutionHistory)
Scalar computeLinearRegressionLogLog(std::vector< Scalar > &x, std::vector< Scalar > &y)
Scalar computeLinearRegression(std::vector< Scalar > &x, std::vector< Scalar > &y)
Teuchos::RCP< LinearRegression< Scalar > > linearRegression()
bool approxZero(Scalar value, Scalar tol=Teuchos::ScalarTraits< Scalar >::sfmin())
Test if value is approximately zero within tolerance.