Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
01_Utilize_Thyra.cpp
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#include <iomanip>
10#include <iostream>
11#include <stdlib.h>
12#include <math.h>
13#include "Teuchos_StandardCatchMacros.hpp"
14
15#include "Thyra_VectorStdOps.hpp"
16#include "Thyra_DefaultSpmdVectorSpace.hpp"
17#include "Thyra_DetachedVectorView.hpp"
18
19
20using namespace std;
21using Teuchos::RCP;
22
205int main(int argc, char *argv[])
206{
207 bool verbose = true;
208 bool success = false;
209 try {
210 // Solution and its time-derivative.
211 int vectorLength = 2; // number state unknowns
212 RCP<const Thyra::VectorSpaceBase<double> > xSpace =
213 Thyra::defaultSpmdVectorSpace<double>(vectorLength);
214
215 RCP<Thyra::VectorBase<double> > x_n = Thyra::createMember(xSpace);
216 RCP<Thyra::VectorBase<double> > xDot_n = Thyra::createMember(xSpace);
217
218 // Initial Conditions
219 double time = 0.0;
220 double epsilon = 1.0e-1;
221 { // Scope to delete DetachedVectorViews
222 Thyra::DetachedVectorView<double> x_n_view(*x_n);
223 x_n_view[0] = 2.0;
224 x_n_view[1] = 0.0;
225 Thyra::DetachedVectorView<double> xDot_n_view(*xDot_n);
226 xDot_n_view[0] = 0.0;
227 xDot_n_view[1] = -2.0/epsilon;
228 }
229
230 // Timestep size
231 double finalTime = 2.0;
232 int nTimeSteps = 2001;
233 const double constDT = finalTime/(nTimeSteps-1);
234
235 // Advance the solution to the next timestep.
236 int n = 0;
237 bool passing = true;
238 cout << n << " " << time << " " << get_ele(*(x_n), 0)
239 << " " << get_ele(*(x_n), 1) << endl;
240 while (passing && time < finalTime && n < nTimeSteps) {
241
242 // Initialize next time step
243 RCP<Thyra::VectorBase<double> > x_np1 = x_n->clone_v(); // at time index n+1
244
245 // Set the timestep and time.
246 double dt = constDT;
247 time = (n+1)*dt;
248
249 // Righthand side evaluation and time-derivative at n.
250 {
251 Thyra::ConstDetachedVectorView<double> x_n_view(*x_n);
252 Thyra::DetachedVectorView<double> xDot_n_view(*xDot_n);
253 xDot_n_view[0] = x_n_view[1];
254 xDot_n_view[1] =
255 ((1.0-x_n_view[0]*x_n_view[0])*x_n_view[1]-x_n_view[0])/epsilon;
256 }
257
258 // Take the timestep - Forward Euler
259 Thyra::V_VpStV(x_np1.ptr(), *x_n, dt, *xDot_n);
260
261 // Test if solution is passing.
262 if ( std::isnan(Thyra::norm(*x_np1)) ) {
263 passing = false;
264 } else {
265 // Promote to next step (n <- n+1).
266 n++;
267 Thyra::V_V(x_n.ptr(), *x_np1);
268 }
269
270 // Output
271 if ( n%100 == 0 )
272 cout << n << " " << time << " " << get_ele(*(x_n), 0)
273 << " " << get_ele(*(x_n), 1) << endl;
274 }
275
276 // Test for regression.
277 RCP<Thyra::VectorBase<double> > x_regress = x_n->clone_v();
278 {
279 Thyra::DetachedVectorView<double> x_regress_view(*x_regress);
280 x_regress_view[0] = -1.59496108218721311;
281 x_regress_view[1] = 0.96359412806611255;
282 }
283
284 RCP<Thyra::VectorBase<double> > x_error = x_n->clone_v();
285 Thyra::V_VmV(x_error.ptr(), *x_n, *x_regress);
286 double x_L2norm_error = Thyra::norm_2(*x_error );
287 double x_L2norm_regress = Thyra::norm_2(*x_regress);
288
289 cout << "Relative L2 Norm of the error (regression) = "
290 << x_L2norm_error/x_L2norm_regress << endl;
291 if ( x_L2norm_error > 1.0e-08*x_L2norm_regress) {
292 passing = false;
293 cout << "FAILED regression constraint!" << endl;
294 }
295
296 RCP<Thyra::VectorBase<double> > x_best = x_n->clone_v();
297 {
298 Thyra::DetachedVectorView<double> x_best_view(*x_best);
299 x_best_view[0] = -1.59496108218721311;
300 x_best_view[1] = 0.96359412806611255;
301 }
302
303 Thyra::V_VmV(x_error.ptr(), *x_n, *x_best);
304 x_L2norm_error = Thyra::norm_2(*x_error);
305 double x_L2norm_best = Thyra::norm_2(*x_best );
306
307 cout << "Relative L2 Norm of the error (best) = "
308 << x_L2norm_error/x_L2norm_best << endl;
309 if ( x_L2norm_error > 0.02*x_L2norm_best) {
310 passing = false;
311 cout << "FAILED best constraint!" << endl;
312 }
313 if (passing) success = true;
314 }
315 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
316
317 if(success)
318 cout << "\nEnd Result: Test Passed!" << std::endl;
319
320 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
321}
int main(int argc, char *argv[])
Description: