Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
02_Use_ModelEvaluator.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
20
21
22using namespace std;
23using Teuchos::RCP;
24
198int main(int argc, char *argv[])
199{
200 bool verbose = true;
201 bool success = false;
202 try {
203 // Construct ModelEvaluator
204 Teuchos::RCP<const Thyra::ModelEvaluator<double> >
205 model = Teuchos::rcp(new VanDerPol_ModelEvaluator_02<double>());
206
207 // Get initial conditions from ModelEvaluator::getNominalValues.
208 RCP<Thyra::VectorBase<double> > x_n =
209 model->getNominalValues().get_x()->clone_v();
210 RCP<Thyra::VectorBase<double> > xDot_n =
211 model->getNominalValues().get_x_dot()->clone_v();
212
213 // Timestep size
214 double time = 0.0;
215 double finalTime = 2.0;
216 int nTimeSteps = 2001;
217 const double constDT = finalTime/(nTimeSteps-1);
218
219 // Advance the solution to the next timestep.
220 int n = 0;
221 bool passing = true;
222 cout << n << " " << time << " " << get_ele(*(x_n), 0)
223 << " " << get_ele(*(x_n), 1) << endl;
224 while (passing && time < finalTime && n < nTimeSteps) {
225
226 // Initialize next time step
227 RCP<Thyra::VectorBase<double> > x_np1 = x_n->clone_v(); // at time index n+1
228
229 // Set the timestep and time.
230 double dt = constDT;
231 time = (n+1)*dt;
232
233 // For explicit ODE formulation, xDot = f(x, t),
234 // xDot is part of the outArgs.
235 auto inArgs = model->createInArgs();
236 auto outArgs = model->createOutArgs();
237 inArgs.set_t(time);
238 inArgs.set_x(x_n);
239 inArgs.set_x_dot(Teuchos::null);
240 outArgs.set_f(xDot_n);
241
242 // Righthand side evaluation and time-derivative at n.
243 model->evalModel(inArgs, outArgs);
244
245 // Take the timestep - Forward Euler
246 Thyra::V_VpStV(x_np1.ptr(), *x_n, dt, *xDot_n);
247
248 // Test if solution is passing.
249 if ( std::isnan(Thyra::norm(*x_np1)) ) {
250 passing = false;
251 } else {
252 // Promote to next step (n <- n+1).
253 n++;
254 Thyra::V_V(x_n.ptr(), *x_np1);
255 }
256
257 // Output
258 if ( n%100 == 0 )
259 cout << n << " " << time << " " << get_ele(*(x_n), 0)
260 << " " << get_ele(*(x_n), 1) << endl;
261 }
262
263 // Test for regression.
264 RCP<Thyra::VectorBase<double> > x_regress = x_n->clone_v();
265 {
266 Thyra::DetachedVectorView<double> x_regress_view(*x_regress);
267 x_regress_view[0] = -1.59496108218721311;
268 x_regress_view[1] = 0.96359412806611255;
269 }
270
271 RCP<Thyra::VectorBase<double> > x_error = x_n->clone_v();
272 Thyra::V_VmV(x_error.ptr(), *x_n, *x_regress);
273 double x_L2norm_error = Thyra::norm_2(*x_error );
274 double x_L2norm_regress = Thyra::norm_2(*x_regress);
275
276 cout << "Relative L2 Norm of the error (regression) = "
277 << x_L2norm_error/x_L2norm_regress << endl;
278 if ( x_L2norm_error > 1.0e-08*x_L2norm_regress) {
279 passing = false;
280 cout << "FAILED regression constraint!" << endl;
281 }
282
283 RCP<Thyra::VectorBase<double> > x_best = x_n->clone_v();
284 {
285 Thyra::DetachedVectorView<double> x_best_view(*x_best);
286 x_best_view[0] = -1.59496108218721311;
287 x_best_view[1] = 0.96359412806611255;
288 }
289
290 Thyra::V_VmV(x_error.ptr(), *x_n, *x_best);
291 x_L2norm_error = Thyra::norm_2(*x_error);
292 double x_L2norm_best = Thyra::norm_2(*x_best );
293
294 cout << "Relative L2 Norm of the error (best) = "
295 << x_L2norm_error/x_L2norm_best << endl;
296 if ( x_L2norm_error > 0.02*x_L2norm_best) {
297 passing = false;
298 cout << "FAILED best constraint!" << endl;
299 }
300 if (passing) success = true;
301 }
302 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
303
304 if(success)
305 cout << "\nEnd Result: Test Passed!" << std::endl;
306
307 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
308}
int main(int argc, char *argv[])
Description: