Tempus  Version of the Day Time Integration
00_Basic_Problem.cpp
Go to the documentation of this file.
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
6 // ****************************************************************************
8
9 #include <iomanip>
10 #include <iostream>
11 #include <stdlib.h>
12 #include <math.h>
13 #include "Teuchos_StandardCatchMacros.hpp"
14
15 using namespace std;
16
216 int 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  int n = 0;
227  double time = 0.0;
228  double epsilon = 1.0e-1;
229  bool passed = true; // ICs are considered passed.
230  x_n [0] = 2.0;
231  x_n [1] = 0.0;
232  xDot_n[0] = 0.0;
233  xDot_n[1] = -2.0/epsilon;
234
235  // Timestep size
236  double finalTime = 2.0;
237  int nTimeSteps = 2001;
238  const double constDT = finalTime/(nTimeSteps-1);
239
240  // Advance the solution to the next timestep.
241  cout << n << " " << time << " " << x_n[0] << " " << x_n[1] << endl;
242  while (passed && 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 passed.
262  if ( std::isnan(x_n[0]) || std::isnan(x_n[1]) ) {
263  passed = false;
264  } else {
265  // Promote to next step (n <- n+1).
266  x_n[0] = x_np1[0];
267  x_n[1] = x_np1[1];
268  n++;
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  passed = 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  passed = false;
311  cout << "FAILED best constraint!" << endl;
312  }
313  if (passed) 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[])
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)