Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
00_Basic_Problem.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
15using namespace std;
16
216int main(int argc, char *argv[])
217{
218 bool verbose = true;
219 bool success = false;
220 try {
221 // Solution and its time-derivative.
222 double x_n[2]; // at time index n
223 double xDot_n[2]; // at time index n
224
225 // Initial Conditions
226 double time = 0.0;
227 double epsilon = 1.0e-1;
228 x_n [0] = 2.0;
229 x_n [1] = 0.0;
230 xDot_n[0] = 0.0;
231 xDot_n[1] = -2.0/epsilon;
232
233 // Timestep size
234 double finalTime = 2.0;
235 int nTimeSteps = 2001;
236 const double constDT = finalTime/(nTimeSteps-1);
237
238 // Advance the solution to the next timestep.
239 int n = 0;
240 bool passing = true;
241 cout << n << " " << time << " " << x_n[0] << " " << x_n[1] << endl;
242 while (passing && time < finalTime && n < nTimeSteps) {
243
244 // Initialize next time step
245 double x_np1[2]; // at time index n+1
246 x_np1[0] = x_n[0];
247 x_np1[1] = x_n[1];
248
249 // Set the timestep and time.
250 double dt = constDT;
251 time = (n+1)*dt;
252
253 // Righthand side evaluation and time-derivative at n.
254 xDot_n[0] = x_n[1];
255 xDot_n[1] = ((1.0 - x_n[0]*x_n[0])*x_n[1] - x_n[0])/epsilon;
256
257 // Take the timestep - Forward Euler
258 x_np1[0] = x_n[0] + dt*xDot_n[0];
259 x_np1[1] = x_n[1] + dt*xDot_n[1];
260
261 // Test if solution is passing.
262 if ( std::isnan(x_n[0]) || std::isnan(x_n[1]) ) {
263 passing = false;
264 } else {
265 // Promote to next step (n <- n+1).
266 n++;
267 x_n[0] = x_np1[0];
268 x_n[1] = x_np1[1];
269 }
270
271 // Output
272 if ( n%100 == 0 )
273 cout << n << " " << time << " " << x_n[0] << " " << x_n[1] << endl;
274
275 }
276
277 // Test for regression.
278 double x_regress[2]; // Regression results for nTimeSteps = 2001
279 x_regress[0] = -1.59496108218721311;
280 x_regress[1] = 0.96359412806611255;
281 double x_L2norm_error = 0.0;
282 double x_L2norm_regress = 0.0;
283 for (int i=0; i < 2; i++) {
284 x_L2norm_error += (x_n[i]-x_regress[i])*(x_n[i]-x_regress[i]);
285 x_L2norm_regress += x_regress[1]*x_regress[1];
286 }
287 x_L2norm_error = sqrt(x_L2norm_error );
288 x_L2norm_regress = sqrt(x_L2norm_regress);
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 double x_best[2]; // Best resolution with nTimeSteps = 2000000001
297 x_best[0] = -1.58184083624543947;
298 x_best[1] = 0.97844890081968072;
299 x_L2norm_error = 0.0;
300 double x_L2norm_best = 0.0;
301 for (int i=0; i < 2; i++) {
302 x_L2norm_error += (x_n[i]-x_best[i])*(x_n[i]-x_best[i]);
303 x_L2norm_best += x_best[1]*x_best[1];
304 }
305 x_L2norm_error = sqrt(x_L2norm_error);
306 x_L2norm_best = sqrt(x_L2norm_best );
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: