9 #ifndef TEMPUS_CONVERGENCE_TEST_UTILS_HPP
10 #define TEMPUS_CONVERGENCE_TEST_UTILS_HPP
18 #include "Tempus_Stepper.hpp"
20 namespace Tempus_Test {
25 template <
class Scalar>
29 void setData(std::vector<Scalar>& x, std::vector<Scalar>& y);
39 std::vector<Scalar>
x_;
40 std::vector<Scalar>
y_;
48 template <
class Scalar>
51 isInitialized_ =
false;
54 template <
class Scalar>
56 std::vector<Scalar>& y)
58 validateXYData_(x, y);
61 isInitialized_ =
true;
65 template <
class Scalar>
67 std::vector<Scalar>& y)
70 x.size() != y.size(), std::logic_error,
71 "x and y data are note the same size for linear regression\n");
73 "Not enough data points for linear regression!\n");
74 int N = Teuchos::as<int>(x.size());
78 for (
int i = 1; i < N; ++i) {
86 template <
class Scalar>
93 template <
class Scalar>
100 template <
class Scalar>
106 int N = Teuchos::as<int>(x_.size());
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];
114 sum1 *= Scalar(-2 * ST::one());
115 sum2 *= Scalar(-2 * ST::one());
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];
125 sum3 *= Scalar(2 * ST::one() / Scalar(N));
126 sum4 *= Scalar(2 * ST::one() / Scalar(N));
128 slope_ = (sum3 + sum1) / (sum4 + sum2);
130 yIntercept_ = ST::zero();
131 for (
int i = 0; i < N; ++i) {
132 yIntercept_ += y_[i] - slope_ * x_[i];
134 yIntercept_ *= Scalar(ST::one() / Scalar(N));
138 template <
class Scalar>
146 template <
class Scalar>
148 Scalar& slope, Scalar& yIntercept)
157 template <
class Scalar>
159 std::vector<Scalar>& y)
162 int N = Teuchos::as<int>(x.size());
163 std::vector<Scalar> xlog;
164 std::vector<Scalar> ylog;
166 for (
int i = 0; i < N; ++i) {
168 xlog.push_back(log(x[i]));
169 ylog.push_back(log(y[i]));
178 template <
class Scalar>
186 template <
class Scalar>
189 std::vector<Scalar>& StepSize,
191 std::vector<Scalar>& xErrorNorm, Scalar& xSlope,
193 std::vector<Scalar>& xDotErrorNorm, Scalar& xDotSlope,
195 std::vector<Scalar>& xDotDotErrorNorm, Scalar& xDotDotSlope,
198 if (solutions.empty()) {
199 out <<
"No solutions to write out!\n"
203 std::size_t lastElem = solutions.size() - 1;
204 std::vector<Scalar> dt;
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);
214 xSlope = computeLinearRegressionLogLog<Scalar>(dt, xErrorNorm);
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);
224 xDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotErrorNorm);
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);
235 xDotDotSlope = computeLinearRegressionLogLog<Scalar>(dt, xDotDotErrorNorm);
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;
249 std::ofstream ftmp(filename);
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];
264 template <
class Scalar>
267 std::vector<Scalar>& StepSize,
269 std::vector<Scalar>& xErrorNorm, Scalar& xSlope,
271 std::vector<Scalar>& xDotErrorNorm, Scalar& xDotSlope,
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);
282 template <
class Scalar>
285 std::vector<Scalar>& StepSize,
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);
300 template <
class Scalar>
302 const std::string filename,
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();
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);
329 template <
class Scalar>
331 const std::string filename,
340 #endif // TEMPUS_CONVERGENCE_TEST_UTILS_HPP
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Scalar getYIntercept() const
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)