Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEulerInt.cpp
Go to the documentation of this file.
1 /**
2  * @file BEulerInt.cpp
3  */
4 
5 /*
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  * See file License.txt for licensing information.
10  */
11 
14 #include "cantera/base/global.h"
15 
16 using namespace std;
17 
18 namespace Cantera
19 {
20 
21 BEulerErr::BEulerErr(const std::string& msg) :
22  CanteraError("BEulerInt", msg)
23 {
24 }
25 
27  m_iter(Newton_Iter),
28  m_method(BEulerVarStep),
29  m_jacFormMethod(BEULER_JAC_NUM),
30  m_rowScaling(true),
31  m_colScaling(false),
32  m_matrixConditioning(false),
33  m_itol(0),
34  m_reltol(1.e-4),
35  m_abstols(1.e-10),
36  m_hmax(0.0),
37  m_maxord(0),
38  m_time_step_num(0),
39  m_time_step_attempts(0),
40  m_max_time_step_attempts(11000000),
41  m_numInitialConstantDeltaTSteps(0),
42  m_failure_counter(0),
43  m_min_newt_its(0),
44  m_printSolnStepInterval(1),
45  m_printSolnNumberToTout(1),
46  m_printSolnFirstSteps(0),
47  m_dumpJacobians(false),
48  m_neq(0),
49  m_t0(0.0),
50  m_time_final(0.0),
51  time_n(0.0),
52  time_nm1(0.0),
53  time_nm2(0.0),
54  delta_t_n(0.0),
55  delta_t_nm1(0.0),
56  delta_t_nm2(0.0),
57  delta_t_np1(1.0E-8),
58  delta_t_max(1.0E300),
59  m_func(0),
60  tdjac_ptr(0),
61  m_print_flag(3),
62  m_nfe(0),
63  m_nJacEval(0),
64  m_numTotalNewtIts(0),
65  m_numTotalLinearSolves(0),
66  m_numTotalConvFails(0),
67  m_numTotalTruncFails(0),
68  num_failures(0)
69 {
70  warn_deprecated("class BEulerInt", "To be removed after Cantera 2.2.");
71 }
72 
73 BEulerInt::~BEulerInt()
74 {
75  delete tdjac_ptr;
76 }
77 
78 void BEulerInt::setTolerances(double reltol, size_t n, double* abstol)
79 {
80  m_itol = 1;
81  m_abstol.resize(m_neq);
82  if (static_cast<int>(n) != m_neq) {
83  printf("ERROR n is wrong\n");
84  exit(-1);
85  }
86  for (int i = 0; i < m_neq; i++) {
87  m_abstol[i] = abstol[i];
88  }
89  m_reltol = reltol;
90 }
91 
92 void BEulerInt::setTolerances(double reltol, double abstol)
93 {
94  m_itol = 0;
95  m_reltol = reltol;
96  m_abstols = abstol;
97 }
98 
99 void BEulerInt::setProblemType(int jacFormMethod)
100 {
101  m_jacFormMethod = jacFormMethod;
102 }
103 
104 void BEulerInt::setMethodBEMT(BEulerMethodType t)
105 {
106  m_method = t;
107 }
108 
109 void BEulerInt::setMaxStep(doublereal hmax)
110 {
111  m_hmax = hmax;
112 }
113 
114 void BEulerInt::setMaxNumTimeSteps(int maxNumTimeSteps)
115 {
116  m_max_time_step_attempts = maxNumTimeSteps;
117 }
118 
119 void BEulerInt::setNumInitialConstantDeltaTSteps(int num)
120 {
122 }
123 
124 void BEulerInt::setPrintSolnOptions(int printSolnStepInterval,
125  int printSolnNumberToTout,
126  int printSolnFirstSteps,
127  bool dumpJacobians)
128 {
129  m_printSolnStepInterval = printSolnStepInterval;
130  m_printSolnNumberToTout = printSolnNumberToTout;
131  m_printSolnFirstSteps = printSolnFirstSteps;
132  m_dumpJacobians = dumpJacobians;
133 }
134 
136 {
137  m_iter = t;
138 }
139 
140 void BEulerInt::setNonLinOptions(int min_newt_its, bool matrixConditioning,
141  bool colScaling, bool rowScaling)
142 {
143  m_min_newt_its = min_newt_its;
144  m_matrixConditioning = matrixConditioning;
145  m_colScaling = colScaling;
146  m_rowScaling = rowScaling;
147  if (m_colScaling && m_colScales.empty()) {
148  m_colScales.assign(m_neq, 1.0);
149  }
150  if (m_rowScaling && m_rowScales.empty()) {
151  m_rowScales.assign(m_neq, 1.0);
152  }
153 }
154 
155 void BEulerInt::setInitialTimeStep(double deltaT)
156 {
157  delta_t_np1 = deltaT;
158 }
159 
160 void BEulerInt::setPrintFlag(int print_flag)
161 {
162  m_print_flag = print_flag;
163 }
164 
166 {
167  m_neq = func.nEquations();
168  m_t0 = t0;
169  internalMalloc();
170 
171  /*
172  * Get the initial conditions.
173  */
174  func.getInitialConditions(m_t0, &m_y_n[0], &m_ydot_n[0]);
175 
176  // Store a pointer to the residual routine in the object
177  m_func = &func;
178 
179  /*
180  * Initialize the various time counters in the object
181  */
182  time_n = t0;
183  time_nm1 = time_n;
184  time_nm2 = time_nm1;
185  delta_t_n = 0.0;
186  delta_t_nm1 = 0.0;
187 }
188 
189 void BEulerInt::reinitializeRJE(double t0, ResidJacEval& func)
190 {
191  m_neq = func.nEquations();
192  m_t0 = t0;
193  internalMalloc();
194  /*
195  * At the initial time, get the initial conditions and time and store
196  * them into internal storage in the object, my[].
197  */
198  m_t0 = t0;
199  func.getInitialConditions(m_t0, &m_y_n[0], &m_ydot_n[0]);
200  /**
201  * Set up the internal weights that are used for testing convergence
202  */
203  setSolnWeights();
204 
205  // Store a pointer to the function
206  m_func = &func;
207 
208 }
209 
210 double BEulerInt::getPrintTime(double time_current)
211 {
212  double tnext;
213  if (m_printSolnNumberToTout > 0) {
214  double dt = (m_time_final - m_t0) / m_printSolnNumberToTout;
215  for (int i = 0; i <= m_printSolnNumberToTout; i++) {
216  tnext = m_t0 + dt * i;
217  if (tnext >= time_current) {
218  return tnext;
219  }
220  }
221  }
222  return 1.0E300;
223 }
224 
225 int BEulerInt::nEvals() const
226 {
227  return m_nfe;
228 }
229 
231 {
232  m_ewt.assign(m_neq, 0.0);
233  m_y_n.assign(m_neq, 0.0);
234  m_y_nm1.assign(m_neq, 0.0);
235  m_y_pred_n.assign(m_neq, 0.0);
236  m_ydot_n.assign(m_neq, 0.0);
237  m_ydot_nm1.assign(m_neq, 0.0);
238  m_resid.assign(m_neq, 0.0);
239  m_residWts.assign(m_neq, 0.0);
240  m_wksp.assign(m_neq, 0.0);
241  if (m_rowScaling) {
242  m_rowScales.assign(m_neq, 1.0);
243  }
244  if (m_colScaling) {
245  m_colScales.assign(m_neq, 1.0);
246  }
248 }
249 
251 {
252  int i;
253  if (m_itol == 1) {
254  for (i = 0; i < m_neq; i++) {
255  m_ewt[i] = m_abstol[i] + m_reltol * 0.5 *
256  (fabs(m_y_n[i]) + fabs(m_y_pred_n[i]));
257  }
258  } else {
259  for (i = 0; i < m_neq; i++) {
260  m_ewt[i] = m_abstols + m_reltol * 0.5 *
261  (fabs(m_y_n[i]) + fabs(m_y_pred_n[i]));
262  }
263  }
264 }
265 
267 {
268  m_func->calcSolnScales(time_n, &m_y_n[0], &m_y_nm1[0], &m_colScales[0]);
269 }
270 
272 {
273  /*
274  * We compute residual weights here, which we define as the L_0 norm
275  * of the Jacobian Matrix, weighted by the solution weights.
276  * This is the proper way to guage the magnitude of residuals. However,
277  * it does need the evaluation of the Jacobian, and the implementation
278  * below is slow, but doesn't take up much memory.
279  *
280  * Here a small weighting indicates that the change in solution is
281  * very sensitive to that equation.
282  */
283  int i, j;
284  double* data = &(*(jac.begin()));
285  double value;
286  for (i = 0; i < m_neq; i++) {
287  m_residWts[i] = fabs(data[i] * m_ewt[0]);
288  for (j = 1; j < m_neq; j++) {
289  value = fabs(data[j*m_neq + i] * m_ewt[j]);
290  m_residWts[i] = std::max(m_residWts[i], value);
291  }
292  }
293 }
294 
295 double BEulerInt::filterNewStep(double timeCurrent, double* y_current, double* ydot_current)
296 {
297  return 0.0;
298 }
299 
300 /*
301  * Print out for relevant time step information
302  */
303 static void print_time_step1(int order, int n_time_step, double time,
304  double delta_t_n, double delta_t_nm1,
305  bool step_failed, int num_failures)
306 {
307  const char* string = 0;
308  if (order == 0) {
309  string = "Backward Euler";
310  } else if (order == 1) {
311  string = "Forward/Backward Euler";
312  } else if (order == 2) {
313  string = "Adams-Bashforth/TR";
314  }
315  writeline('=', 80, true, true);
316  printf("\nStart of Time Step: %5d Time_n = %9.5g Time_nm1 = %9.5g\n",
317  n_time_step, time, time - delta_t_n);
318  printf("\tIntegration method = %s\n", string);
319  if (step_failed) {
320  printf("\tPreviously attempted step was a failure\n");
321  }
322  if (delta_t_n > delta_t_nm1) {
323  string = "(Increased from previous iteration)";
324  } else if (delta_t_n < delta_t_nm1) {
325  string = "(Decreased from previous iteration)";
326  } else {
327  string = "(same as previous iteration)";
328  }
329  printf("\tdelta_t_n = %8.5e %s", delta_t_n, string);
330  if (num_failures > 0) {
331  printf("\t(Bad_History Failure Counter = %d)", num_failures);
332  }
333  printf("\n\tdelta_t_nm1 = %8.5e\n", delta_t_nm1);
334 }
335 
336 /*
337  * Print out for relevant time step information
338  */
339 static void print_time_step2(int time_step_num, int order,
340  double time, double time_error_factor,
341  double delta_t_n, double delta_t_np1)
342 {
343  printf("\tTime Step Number %5d was a success: time = %10g\n", time_step_num,
344  time);
345  printf("\t\tEstimated Error\n");
346  printf("\t\t-------------------- = %8.5e\n", time_error_factor);
347  printf("\t\tTolerated Error\n\n");
348  printf("\t- Recommended next delta_t (not counting history) = %g\n",
349  delta_t_np1);
350  writeline('=', 80, true, true);
351 }
352 
353 /*
354  * Print Out descriptive information on why the current step failed
355  */
356 static void print_time_fail(bool convFailure, int time_step_num,
357  double time, double delta_t_n,
358  double delta_t_np1, double time_error_factor)
359 {
360  writeline('=', 80, true, true);
361  if (convFailure) {
362  printf("\tTime Step Number %5d experienced a convergence "
363  "failure\n", time_step_num);
364  printf("\tin the non-linear or linear solver\n");
365  printf("\t\tValue of time at failed step = %g\n", time);
366  printf("\t\tdelta_t of the failed step = %g\n",
367  delta_t_n);
368  printf("\t\tSuggested value of delta_t to try next = %g\n",
369  delta_t_np1);
370  } else {
371  printf("\tTime Step Number %5d experienced a truncation error "
372  "failure!\n", time_step_num);
373  printf("\t\tValue of time at failed step = %g\n", time);
374  printf("\t\tdelta_t of the failed step = %g\n",
375  delta_t_n);
376  printf("\t\tSuggested value of delta_t to try next = %g\n",
377  delta_t_np1);
378  printf("\t\tCalculated truncation error factor = %g\n",
379  time_error_factor);
380  }
381  writeline('=', 80, true, true);
382 }
383 
384 /*
385  * Print out the final results and counters
386  */
387 static void print_final(double time, int step_failed,
388  int time_step_num, int num_newt_its,
389  int total_linear_solves, int numConvFails,
390  int numTruncFails, int nfe, int nJacEval)
391 {
392  writeline('=', 80, true, true);
393  printf("TIME INTEGRATION ROUTINE HAS FINISHED: ");
394  if (step_failed) {
395  printf(" IT WAS A FAILURE\n");
396  } else {
397  printf(" IT WAS A SUCCESS\n");
398  }
399  printf("\tEnding time = %g\n", time);
400  printf("\tNumber of time steps = %d\n", time_step_num);
401  printf("\tNumber of newt its = %d\n", num_newt_its);
402  printf("\tNumber of linear solves = %d\n", total_linear_solves);
403  printf("\tNumber of convergence failures= %d\n", numConvFails);
404  printf("\tNumber of TimeTruncErr fails = %d\n", numTruncFails);
405  printf("\tNumber of Function evals = %d\n", nfe);
406  printf("\tNumber of Jacobian evals/solvs= %d\n", nJacEval);
407  writeline('=', 80, true, true);
408 }
409 
410 /*
411  * Header info for one line comment about a time step
412  */
413 static void print_lvl1_Header(int nTimes)
414 {
415  printf("\n");
416  if (nTimes) {
417  writeline('-', 80);
418  }
419  printf("time Time Time Time ");
420  if (nTimes == 0) {
421  printf(" START");
422  } else {
423  printf(" (continued)");
424  }
425  printf("\n");
426 
427  printf("step (sec) step Newt Aztc bktr trunc ");
428  printf("\n");
429 
430  printf(" No. Rslt size Its Its stps error |");
431  printf(" comment");
432  writeline('-', 80, true);
433 }
434 
435 /*
436  * One line entry about time step
437  * rslt -> 4 letter code
438  */
439 static void print_lvl1_summary(
440  int time_step_num, double time, const char* rslt, double delta_t_n,
441  int newt_its, int aztec_its, int bktr_stps, double time_error_factor,
442  const char* comment)
443 {
444  printf("%6d %11.6g %4s %10.4g %4d %4d %4d %11.4g",
445  time_step_num, time, rslt, delta_t_n, newt_its, aztec_its,
446  bktr_stps, time_error_factor);
447  if (comment) {
448  printf(" | %s", comment);
449  }
450  printf("\n");
451 }
452 
453 /*
454  * This routine subtracts 2 numbers. If the difference is less
455  * than 1.0E-14 times the magnitude of the smallest number,
456  * then diff returns an exact zero.
457  * It also returns an exact zero if the difference is less than
458  * 1.0E-300.
459  *
460  * returns: a - b
461  *
462  * This routine is used in numerical differencing schemes in order
463  * to avoid roundoff errors resulting in creating Jacobian terms.
464  * Note: This is a slow routine. However, Jacobian errors may cause
465  * loss of convergence. Therefore, in practice this routine
466  * has proved cost-effective.
467  */
468 double subtractRD(double a, double b)
469 {
470  double diff = a - b;
471  double d = std::min(fabs(a), fabs(b));
472  d *= 1.0E-14;
473  double ad = fabs(diff);
474  if (ad < 1.0E-300) {
475  diff = 0.0;
476  }
477  if (ad < d) {
478  diff = 0.0;
479  }
480  return diff;
481 }
482 
483 void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f,
484  double time_curr, double CJ,
485  double* const y,
486  double* const ydot,
487  int num_newt_its)
488 {
489  int i, j;
490  double* col_j;
491  double ysave, ydotsave, dy;
492  /**
493  * Clear the factor flag
494  */
495  J.clearFactorFlag();
496 
497 
498  if (m_jacFormMethod & BEULER_JAC_ANAL) {
499  /********************************************************************
500  * Call the function to get a Jacobian.
501  */
502  m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
503  m_nJacEval++;
504  m_nfe++;
505  } else {
506  /*******************************************************************
507  * Generic algorithm to calculate a numerical Jacobian
508  */
509  /*
510  * Calculate the current value of the RHS given the
511  * current conditions.
512  */
513 
514  m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
515  m_nfe++;
516  m_nJacEval++;
517 
518 
519  /*
520  * Malloc a vector and call the function object to return a set of
521  * deltaY's that are appropriate for calculating the numerical
522  * derivative.
523  */
524  vector_fp dyVector(m_neq);
525  m_func->calcDeltaSolnVariables(time_curr, y, &m_y_nm1[0], &dyVector[0],
526  &m_ewt[0]);
527 #ifdef DEBUG_HKM
528  bool print_NumJac = false;
529  if (print_NumJac) {
530  FILE* idy = fopen("NumJac.csv", "w");
531  fprintf(idy, "Unk m_ewt y "
532  "dyVector ResN\n");
533  for (int iii = 0; iii < m_neq; iii++) {
534  fprintf(idy, " %4d %16.8e %16.8e %16.8e %16.8e \n",
535  iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
536  }
537  fclose(idy);
538  }
539 #endif
540  /*
541  * Loop over the variables, formulating a numerical derivative
542  * of the dense matrix.
543  * For the delta in the variable, we will use a variety of approaches
544  * The original approach was to use the error tolerance amount.
545  * This may not be the best approach, as it could be overly large in
546  * some instances and overly small in others.
547  * We will first protect from being overly small, by using the usual
548  * sqrt of machine precision approach, i.e., 1.0E-7,
549  * to bound the lower limit of the delta.
550  */
551  for (j = 0; j < m_neq; j++) {
552 
553 
554  /*
555  * Get a pointer into the column of the matrix
556  */
557 
558 
559  col_j = (double*) J.ptrColumn(j);
560  ysave = y[j];
561  dy = dyVector[j];
562 
563  y[j] = ysave + dy;
564  dy = y[j] - ysave;
565  ydotsave = ydot[j];
566  ydot[j] += dy * CJ;
567  /*
568  * Call the functon
569  */
570 
571 
572  m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, &m_wksp[0],
573  JacDelta_ResidEval, j, dy);
574  m_nfe++;
575  double diff;
576  for (i = 0; i < m_neq; i++) {
577  diff = subtractRD(m_wksp[i], f[i]);
578  col_j[i] = diff / dy;
579  }
580 
581  y[j] = ysave;
582  ydot[j] = ydotsave;
583 
584  }
585  }
586 
587 
588 }
589 
590 void BEulerInt::calc_y_pred(int order)
591 {
592  int i;
593  double c1, c2;
594  switch (order) {
595  case 0:
596  case 1:
597  c1 = delta_t_n;
598  for (i = 0; i < m_neq; i++) {
599  m_y_pred_n[i] = m_y_n[i] + c1 * m_ydot_n[i];
600  }
601  break;
602  case 2:
603  c1 = delta_t_n * (2.0 + delta_t_n / delta_t_nm1) / 2.0;
604  c2 = (delta_t_n * delta_t_n) / (delta_t_nm1 * 2.0);
605  for (i = 0; i < m_neq; i++) {
606  m_y_pred_n[i] = m_y_n[i] + c1 * m_ydot_n[i] - c2 * m_ydot_nm1[i];
607  }
608  break;
609  }
610 
611  /*
612  * Filter the predictions.
613  */
614  m_func->filterSolnPrediction(time_n, &m_y_pred_n[0]);
615 
616 }
617 
618 void BEulerInt::calc_ydot(int order, double* y_curr, double* ydot_curr)
619 {
620  int i;
621  double c1;
622  switch (order) {
623  case 0:
624  case 1: /* First order forward Euler/backward Euler */
625  c1 = 1.0 / delta_t_n;
626  for (i = 0; i < m_neq; i++) {
627  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]);
628  }
629  return;
630  case 2: /* Second order Adams-Bashforth / Trapezoidal Rule */
631  c1 = 2.0 / delta_t_n;
632  for (i = 0; i < m_neq; i++) {
633  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i];
634  }
635  return;
636  }
637 }
638 
640 {
641  int i;
642  double rel_norm, error;
643 #ifdef DEBUG_HKM
644 #define NUM_ENTRIES 5
645  if (m_print_flag > 2) {
646  int imax[NUM_ENTRIES], j, jnum;
647  double dmax;
648  bool used;
649  printf("\t\ttime step truncation error contributors:\n");
650  printf("\t\t I entry actual predicted "
651  " weight ydot\n");
652  printf("\t\t");
653  writeline('-', 70);
654  for (j = 0; j < NUM_ENTRIES; j++) {
655  imax[j] = -1;
656  }
657  for (jnum = 0; jnum < NUM_ENTRIES; jnum++) {
658  dmax = -1.0;
659  for (i = 0; i < m_neq; i++) {
660  used = false;
661  for (j = 0; j < jnum; j++) {
662  if (imax[j] == i) {
663  used = true;
664  }
665  }
666  if (!used) {
667  error = (m_y_n[i] - m_y_pred_n[i]) / m_ewt[i];
668  rel_norm = sqrt(error * error);
669  if (rel_norm > dmax) {
670  imax[jnum] = i;
671  dmax = rel_norm;
672  }
673  }
674  }
675  if (imax[jnum] >= 0) {
676  i = imax[jnum];
677  printf("\t\t%4d %12.4e %12.4e %12.4e %12.4e %12.4e\n",
678  i, dmax, m_y_n[i], m_y_pred_n[i], m_ewt[i], m_ydot_n[i]);
679  }
680  }
681  printf("\t\t");
682  writeline('-', 70);
683  }
684 #endif
685  rel_norm = 0.0;
686  for (i = 0; i < m_neq; i++) {
687  error = (m_y_n[i] - m_y_pred_n[i]) / m_ewt[i];
688  rel_norm += (error * error);
689  }
690  return sqrt(rel_norm / m_neq);
691 }
692 
693 double BEulerInt::time_step_control(int order, double time_error_factor)
694 {
695  double factor = 0.0, power = 0.0, delta_t;
696  const char* yo = "time_step_control";
697 
698  /*
699  * Special case time_error_factor so that zeroes don't cause a problem.
700  */
701  time_error_factor = std::max(1.0E-50, time_error_factor);
702 
703  /*
704  * Calculate the factor for the change in magnitude of time step.
705  */
706  switch (order) {
707  case 1:
708  factor = 1.0/(2.0 *(time_error_factor));
709  power = 0.5;
710  break;
711  case 2:
712  factor = 1.0/(3.0 * (1.0 + delta_t_nm1 / delta_t_n)
713  * (time_error_factor));
714  power = 0.3333333333333333;
715  }
716  factor = pow(factor, power);
717  if (factor < 0.5) {
718  if (m_print_flag > 1) {
719  printf("\t%s: WARNING - Current time step will be chucked\n", yo);
720  printf("\t\tdue to a time step truncation error failure.\n");
721  }
722  delta_t = - 0.5 * delta_t_n;
723  } else {
724  factor = std::min(factor, 1.5);
725  delta_t = factor * delta_t_n;
726  }
727  return delta_t;
728 }
729 
730 double BEulerInt::integrateRJE(double tout, double time_init)
731 {
732  double time_current;
733  bool weAreNotFinished = true;
734  m_time_final = tout;
735  int flag = SUCCESS;
736  /**
737  * Initialize the time step number to zero. step will increment so that
738  * the first time step is number 1
739  */
740  m_time_step_num = 0;
741 
742 
743  /*
744  * Do the integration a step at a time
745  */
746  int istep = 0;
747  int printStep = 0;
748  bool doPrintSoln = false;
749  time_current = time_init;
750  time_n = time_init;
751  time_nm1 = time_init;
752  time_nm2 = time_init;
753  m_func->evalTimeTrackingEqns(time_current, 0.0, &m_y_n[0], &m_ydot_n[0]);
754  double print_time = getPrintTime(time_current);
755  if (print_time == time_current) {
756  m_func->writeSolution(4, time_current, delta_t_n,
757  istep, &m_y_n[0], &m_ydot_n[0]);
758  }
759  /*
760  * We print out column headers here for the case of
761  */
762  if (m_print_flag == 1) {
763  print_lvl1_Header(0);
764  }
765  /*
766  * Call a different user routine at the end of each step,
767  * that will probably print to a file.
768  */
769  m_func->user_out2(0, time_current, 0.0, &m_y_n[0], &m_ydot_n[0]);
770 
771  do {
772 
773  print_time = getPrintTime(time_current);
774  if (print_time >= tout) {
775  print_time = tout;
776  }
777 
778  /************************************************************
779  * Step the solution
780  */
781  time_current = step(tout);
782  istep++;
783  printStep++;
784  /***********************************************************/
785  if (time_current < 0.0) {
786  if (time_current == -1234.) {
787  time_current = 0.0;
788  } else {
789  time_current = -time_current;
790  }
791  flag = FAILURE;
792  }
793 
794  if (flag != FAILURE) {
795  bool retn =
796  m_func->evalStoppingCritera(time_current, delta_t_n,
797  &m_y_n[0], &m_ydot_n[0]);
798  if (retn) {
799  weAreNotFinished = false;
800  doPrintSoln = true;
801  }
802  }
803 
804  /*
805  * determine conditional printing of soln
806  */
807  if (time_current >= print_time) {
808  doPrintSoln = true;
809  }
810  if (m_printSolnStepInterval == printStep) {
811  doPrintSoln = true;
812  }
813  if (m_printSolnFirstSteps > istep) {
814  doPrintSoln = true;
815  }
816 
817  /*
818  * Evaluate time integrated quantities that are calculated at the
819  * end of every successful time step.
820  */
821  if (flag != FAILURE) {
822  m_func->evalTimeTrackingEqns(time_current, delta_t_n,
823  &m_y_n[0], &m_ydot_n[0]);
824  }
825 
826  /*
827  * Call the printout routine.
828  */
829  if (doPrintSoln) {
830  m_func->writeSolution(1, time_current, delta_t_n,
831  istep, &m_y_n[0], &m_ydot_n[0]);
832  printStep = 0;
833  doPrintSoln = false;
834  if (m_print_flag == 1) {
835  print_lvl1_Header(1);
836  }
837  }
838  /*
839  * Call a different user routine at the end of each step,
840  * that will probably print to a file.
841  */
842  if (flag == FAILURE) {
843  m_func->user_out2(-1, time_current, delta_t_n, &m_y_n[0], &m_ydot_n[0]);
844  } else {
845  m_func->user_out2(1, time_current, delta_t_n, &m_y_n[0], &m_ydot_n[0]);
846  }
847 
848  } while (time_current < tout &&
849  m_time_step_attempts < m_max_time_step_attempts &&
850  flag == SUCCESS && weAreNotFinished);
851 
852  /*
853  * Check current time against the max solution time.
854  */
855  if (time_current >= tout) {
856  printf("Simulation completed time integration in %d time steps\n",
858  printf("Final Time: %e\n\n", time_current);
859  } else if (m_time_step_attempts >= m_max_time_step_attempts) {
860  printf("Simulation ran into time step attempt limit in"
861  "%d time steps\n",
863  printf("Final Time: %e\n\n", time_current);
864  } else if (flag == FAILURE) {
865  printf("ERROR: time stepper failed at time = %g\n", time_current);
866  }
867 
868  /*
869  * Print out the final results and counters.
870  */
871  print_final(time_n, flag, m_time_step_num, m_numTotalNewtIts,
874 
875  /*
876  * Call a different user routine at the end of each step,
877  * that will probably print to a file.
878  */
879  m_func->user_out2(2, time_current, delta_t_n, &m_y_n[0], &m_ydot_n[0]);
880 
881 
882  if (flag != SUCCESS) {
883  throw BEulerErr(" BEuler error encountered.");
884  }
885  return time_current;
886 }
887 
888 double BEulerInt::step(double t_max)
889 {
890  double CJ;
891  bool step_failed = false;
892  bool giveUp = false;
893  bool convFailure = false;
894  const char* rslt;
895  double time_error_factor = 0.0;
896  double normFilter = 0.0;
897  int numTSFailures = 0;
898  int bktr_stps = 0;
899  int nonlinearloglevel = m_print_flag;
900  int num_newt_its = 0;
901  int aztec_its = 0;
902  string comment;
903  /*
904  * Increment the time counter - May have to be taken back,
905  * if time step is found to be faulty.
906  */
907  m_time_step_num++;
908 
909  /**
910  * Loop here until we achieve a successful step or we set the giveUp
911  * flag indicating that repeated errors have occurred.
912  */
913  do {
914  m_time_step_attempts++;
915  comment.clear();
916 
917  /*
918  * Possibly adjust the delta_t_n value for this time step from the
919  * recommended delta_t_np1 value determined in the previous step
920  * due to maximum time step constraints or other occurences,
921  * known to happen at a given time.
922  */
923  if ((time_n + delta_t_np1) >= t_max) {
924  delta_t_np1 =t_max - time_n;
925  }
926 
927  if (delta_t_np1 >= delta_t_max) {
928  delta_t_np1 = delta_t_max;
929  }
930 
931  /*
932  * Increment the delta_t counters and the time for the current
933  * time step.
934  */
935 
936  delta_t_nm2 = delta_t_nm1;
937  delta_t_nm1 = delta_t_n;
938  delta_t_n = delta_t_np1;
939  time_n += delta_t_n;
940 
941  /*
942  * Determine the integration order of the current step.
943  *
944  * Special case for start-up of time integration procedure
945  * First time step = Do a predictor step as we
946  * have recently added an initial
947  * ydot input option. And, setting ydot=0
948  * is equivalent to not doing a
949  * predictor step.
950  * Second step = If 2nd order method, do a first order
951  * step for this time-step, only.
952  *
953  * If 2nd order method with a constant time step, the
954  * first and second steps are 1/10 the specified step, and
955  * the third step is 8/10 the specified step. This reduces
956  * the error asociated with using lower order
957  * integration on the first two steps. (RCS 11-6-97)
958  *
959  * If the previous time step failed for one reason or another,
960  * do a linear step. It's more robust.
961  */
962  if (m_time_step_num == 1) {
963  m_order = 1; /* Backward Euler */
964  } else if (m_time_step_num == 2) {
965  m_order = 1; /* Forward/Backward Euler */
966  } else if (step_failed) {
967  m_order = 1; /* Forward/Backward Euler */
968  } else if (m_time_step_num > 2) {
969  m_order = 1; /* Specified
970  Predictor/Corrector
971  - not implemented */
972  }
973 
974  /*
975  * Print out an initial statement about the step.
976  */
977  if (m_print_flag > 1) {
978  print_time_step1(m_order, m_time_step_num, time_n, delta_t_n,
979  delta_t_nm1, step_failed, m_failure_counter);
980  }
981 
982  /*
983  * Calculate the predicted solution, m_y_pred_n, for the current
984  * time step.
985  */
987 
988  /*
989  * HKM - Commented this out. I may need it for particles later.
990  * If Solution bounds checking is turned on, we need to crop the
991  * predicted solution to make sure bounds are enforced
992  *
993  *
994  * cropNorm = 0.0;
995  * if (Cur_Realm->Realm_Nonlinear.Constraint_Backtracking_Flag ==
996  * Constraint_Backtrack_Enable) {
997  * cropNorm = cropPredictor(mesh, x_pred_n, abs_time_error,
998  * m_reltol);
999  */
1000 
1001  /*
1002  * Save the old solution, before overwriting with the new solution
1003  * - use
1004  */
1005  m_y_nm1 = m_y_n;
1006 
1007  /*
1008  * Use the predicted value as the initial guess for the corrector
1009  * loop, for
1010  * every step other than the first step.
1011  */
1012  if (m_order > 0) {
1013  m_y_n = m_y_pred_n;
1014  }
1015 
1016  /*
1017  * Save the old time derivative, if necessary, before it is
1018  * overwritten.
1019  * This overwrites ydot_nm1, losing information from the previous time
1020  * step.
1021  */
1022  m_ydot_nm1 = m_ydot_n;
1023 
1024  /*
1025  * Calculate the new time derivative, ydot_n, that is consistent
1026  * with the
1027  * initial guess for the corrected solution vector.
1028  *
1029  */
1030  calc_ydot(m_order, &m_y_n[0], &m_ydot_n[0]);
1031 
1032  /*
1033  * Calculate CJ, the coefficient for the Jacobian corresponding to the
1034  * derivative of the residual wrt to the acceleration vector.
1035  */
1036  if (m_order < 2) {
1037  CJ = 1.0 / delta_t_n;
1038  } else {
1039  CJ = 2.0 / delta_t_n;
1040  }
1041 
1042  /*
1043  * Calculate a new Solution Error Weighting vector
1044  */
1045  setSolnWeights();
1046 
1047  /*
1048  * Solve the system of equations at the current time step.
1049  * Note - x_corr_n and x_dot_n are considered to be updated,
1050  * on return from this solution.
1051  */
1052  int ierror = solve_nonlinear_problem(&m_y_n[0], &m_ydot_n[0],
1053  CJ, time_n, *tdjac_ptr, num_newt_its,
1054  aztec_its, bktr_stps,
1055  nonlinearloglevel);
1056  /*
1057  * Set the appropriate flags if a convergence failure is detected.
1058  */
1059  if (ierror < 0) { /* Step failed */
1060  convFailure = true;
1061  step_failed = true;
1062  rslt = "fail";
1064  m_failure_counter +=3;
1065  if (m_print_flag > 1) {
1066  printf("\tStep is Rejected, nonlinear problem didn't converge,"
1067  "ierror = %d\n", ierror);
1068  }
1069  } else { /* Step succeeded */
1070  convFailure = false;
1071  step_failed = false;
1072  rslt = "done";
1073 
1074  /*
1075  * Apply a filter to a new successful step
1076  */
1077  normFilter = filterNewStep(time_n, &m_y_n[0], &m_ydot_n[0]);
1078  if (normFilter > 1.0) {
1079  convFailure = true;
1080  step_failed = true;
1081  rslt = "filt";
1082  if (m_print_flag > 1) {
1083  printf("\tStep is Rejected, too large filter adjustment = %g\n",
1084  normFilter);
1085  }
1086  } else if (normFilter > 0.0) {
1087  if (normFilter > 0.3) {
1088  if (m_print_flag > 1) {
1089  printf("\tStep was filtered, norm = %g, next "
1090  "time step adjusted\n", normFilter);
1091  }
1092  } else {
1093  if (m_print_flag > 1) {
1094  printf("\tStep was filtered, norm = %g\n", normFilter);
1095  }
1096  }
1097  }
1098  }
1099 
1100  /*
1101  * Calculate the time step truncation error for the current step.
1102  */
1103  if (!step_failed) {
1104  time_error_factor = time_error_norm();
1105  } else {
1106  time_error_factor = 1000.;
1107  }
1108 
1109  /*
1110  * Dynamic time step control- delta_t_n, delta_t_nm1 are set here.
1111  */
1112  if (step_failed) {
1113  /*
1114  * For convergence failures, decrease the step-size by a factor of
1115  * 4 and try again.
1116  */
1117  delta_t_np1 = 0.25 * delta_t_n;
1118  } else if (m_method == BEulerVarStep) {
1119 
1120  /*
1121  * If we are doing a predictor/corrector method, and we are
1122  * past a certain number of time steps given by the input file
1123  * then either correct the DeltaT for the next time step or
1124  *
1125  */
1126  if ((m_order > 0) &&
1128  delta_t_np1 = time_step_control(m_order, time_error_factor);
1129  if (normFilter > 0.1) {
1130  if (delta_t_np1 > delta_t_n) {
1131  delta_t_np1 = delta_t_n;
1132  }
1133  }
1134 
1135  /*
1136  * Check for Current time step failing due to violation of
1137  * time step
1138  * truncation bounds.
1139  */
1140  if (delta_t_np1 < 0.0) {
1142  step_failed = true;
1143  delta_t_np1 = -delta_t_np1;
1144  m_failure_counter += 2;
1145  comment += "TIME TRUNC FAILURE";
1146  rslt = "TRNC";
1147  }
1148 
1149  /*
1150  * Prevent churning of the time step by not increasing the
1151  * time step,
1152  * if the recent "History" of the time step behavior is still bad
1153  */
1154  else if (m_failure_counter > 0) {
1155  delta_t_np1 = std::min(delta_t_np1, delta_t_n);
1156  }
1157  } else {
1158  delta_t_np1 = delta_t_n;
1159  }
1160 
1161  /* Decrease time step if a lot of Newton Iterations are
1162  * taken.
1163  * The idea being if more or less Newton iteration are taken
1164  * than the target number of iterations, then adjust the time
1165  * step downwards so that the target number of iterations or lower
1166  * is achieved. This
1167  * should prevent step failure by too many Newton iterations because
1168  * the time step becomes too large. CCO
1169  * hkm -> put in num_new_its min of 3 because the time step
1170  * was being altered even when num_newt_its == 1
1171  */
1172  int max_Newton_steps = 10000;
1173  int target_num_iter = 5;
1174  if (num_newt_its > 3000 && !step_failed) {
1175  if (max_Newton_steps != target_num_iter) {
1176  double iter_diff = num_newt_its - target_num_iter;
1177  double iter_adjust_zone = max_Newton_steps - target_num_iter;
1178  double target_time_step = delta_t_n
1179  *(1.0 - iter_diff*fabs(iter_diff)/
1180  ((2.0*iter_adjust_zone*iter_adjust_zone)));
1181  target_time_step = std::max(0.5*delta_t_n, target_time_step);
1182  if (target_time_step < delta_t_np1) {
1183  printf("\tNext time step will be decreased from %g to %g"
1184  " because of new its restraint\n",
1185  delta_t_np1, target_time_step);
1186  delta_t_np1 = target_time_step;
1187  }
1188  }
1189  }
1190 
1191 
1192  }
1193 
1194  /*
1195  * The final loop in the time stepping algorithm depends on whether the
1196  * current step was a success or not.
1197  */
1198  if (step_failed) {
1199  /*
1200  * Increment the counter indicating the number of consecutive
1201  * failures
1202  */
1203  numTSFailures++;
1204  /*
1205  * Print out a statement about the failure of the time step.
1206  */
1207  if (m_print_flag > 1) {
1208  print_time_fail(convFailure, m_time_step_num, time_n, delta_t_n,
1209  delta_t_np1, time_error_factor);
1210  } else if (m_print_flag == 1) {
1211  print_lvl1_summary(m_time_step_num, time_n, rslt, delta_t_n,
1212  num_newt_its, aztec_its, bktr_stps,
1213  time_error_factor,
1214  comment.c_str());
1215  }
1216 
1217  /*
1218  * Change time step counters back to the previous step before
1219  * the failed
1220  * time step occurred.
1221  */
1222  time_n -= delta_t_n;
1223  delta_t_n = delta_t_nm1;
1224  delta_t_nm1 = delta_t_nm2;
1225 
1226  /*
1227  * Replace old solution vector and time derivative solution vector.
1228  */
1229  m_y_n = m_y_nm1;
1230  m_ydot_n = m_ydot_nm1;
1231  /*
1232  * Decide whether to bail on the whole loop
1233  */
1234  if (numTSFailures > 35) {
1235  giveUp = true;
1236  }
1237  }
1238 
1239  /*
1240  * Do processing for a successful step.
1241  */
1242  else {
1243 
1244  /*
1245  * Decrement the number of consequative failure counter.
1246  */
1247  m_failure_counter = std::max(0, m_failure_counter-1);
1248 
1249  /*
1250  * Print out final results of a successfull time step.
1251  */
1252  if (m_print_flag > 1) {
1253  print_time_step2(m_time_step_num, m_order, time_n, time_error_factor,
1254  delta_t_n, delta_t_np1);
1255  } else if (m_print_flag == 1) {
1256  print_lvl1_summary(m_time_step_num, time_n, " ", delta_t_n,
1257  num_newt_its, aztec_its, bktr_stps, time_error_factor,
1258  comment.c_str());
1259  }
1260 
1261  /*
1262  * Output information at the end of every successful time step, if
1263  * requested.
1264  *
1265  * fill in
1266  */
1267 
1268 
1269  }
1270  } while (step_failed && !giveUp);
1271 
1272  /*
1273  * Send back the overall result of the time step.
1274  */
1275  if (step_failed) {
1276  if (time_n == 0.0) {
1277  return -1234.0;
1278  }
1279  return -time_n;
1280  }
1281  return time_n;
1282 }
1283 
1284 //-----------------------------------------------------------
1285 // Constants
1286 //-----------------------------------------------------------
1287 
1288 const double DampFactor = 4;
1289 const int NDAMP = 10;
1290 
1291 //-----------------------------------------------------------
1292 // MultiNewton methods
1293 //-----------------------------------------------------------
1294 
1295 double BEulerInt::soln_error_norm(const double* const delta_y,
1296  bool printLargest)
1297 {
1298  int i;
1299  double sum_norm = 0.0, error;
1300  for (i = 0; i < m_neq; i++) {
1301  error = delta_y[i] / m_ewt[i];
1302  sum_norm += (error * error);
1303  }
1304  sum_norm = sqrt(sum_norm / m_neq);
1305  if (printLargest) {
1306  const int num_entries = 8;
1307  double dmax1, normContrib;
1308  int j;
1309  vector_int imax(num_entries, -1);
1310  printf("\t\tPrintout of Largest Contributors to norm "
1311  "of value (%g)\n", sum_norm);
1312  printf("\t\t I ysoln deltaY weightY "
1313  "Error_Norm**2\n");
1314  printf("\t\t ");
1315  writeline('-', 80);
1316  for (int jnum = 0; jnum < num_entries; jnum++) {
1317  dmax1 = -1.0;
1318  for (i = 0; i < m_neq; i++) {
1319  bool used = false;
1320  for (j = 0; j < jnum; j++) {
1321  if (imax[j] == i) {
1322  used = true;
1323  }
1324  }
1325  if (!used) {
1326  error = delta_y[i] / m_ewt[i];
1327  normContrib = sqrt(error * error);
1328  if (normContrib > dmax1) {
1329  imax[jnum] = i;
1330  dmax1 = normContrib;
1331  }
1332  }
1333  }
1334  i = imax[jnum];
1335  if (i >= 0) {
1336  printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e\n",
1337  i, m_y_n[i], delta_y[i], m_ewt[i], dmax1);
1338  }
1339  }
1340  printf("\t\t ");
1341  writeline('-', 80);
1342  }
1343  return sum_norm;
1344 }
1345 #ifdef DEBUG_HKM_JAC
1346 SquareMatrix jacBack();
1347 #endif
1348 
1349 void BEulerInt::doNewtonSolve(double time_curr, double* y_curr,
1350  double* ydot_curr, double* delta_y,
1351  GeneralMatrix& jac, int loglevel)
1352 {
1353  int irow, jcol;
1354 
1355  m_func->evalResidNJ(time_curr, delta_t_n, y_curr,
1356  ydot_curr, delta_y, Base_ResidEval);
1357  m_nfe++;
1358  int sz = m_func->nEquations();
1359  for (int n = 0; n < sz; n++) {
1360  delta_y[n] = -delta_y[n];
1361  }
1362 
1363 
1364  /*
1365  * Column scaling -> We scale the columns of the Jacobian
1366  * by the nominal important change in the solution vector
1367  */
1368  if (m_colScaling) {
1369  if (!jac.factored()) {
1370  /*
1371  * Go get new scales
1372  */
1373  setColumnScales();
1374 
1375  /*
1376  * Scale the new Jacobian
1377  */
1378  double* jptr = &(*(jac.begin()));
1379  for (jcol = 0; jcol < m_neq; jcol++) {
1380  for (irow = 0; irow < m_neq; irow++) {
1381  *jptr *= m_colScales[jcol];
1382  jptr++;
1383  }
1384  }
1385  }
1386  }
1387 
1388  if (m_matrixConditioning) {
1389  if (jac.factored()) {
1390  m_func->matrixConditioning(0, sz, delta_y);
1391  } else {
1392  double* jptr = &(*(jac.begin()));
1393  m_func->matrixConditioning(jptr, sz, delta_y);
1394  }
1395  }
1396 
1397  /*
1398  * row sum scaling -> Note, this is an unequivocal success
1399  * at keeping the small numbers well balanced and
1400  * nonnegative.
1401  */
1402  if (m_rowScaling) {
1403  if (! jac.factored()) {
1404  /*
1405  * Ok, this is ugly. jac.begin() returns an vector<double> iterator
1406  * to the first data location.
1407  * Then &(*()) reverts it to a double *.
1408  */
1409  double* jptr = &(*(jac.begin()));
1410  for (irow = 0; irow < m_neq; irow++) {
1411  m_rowScales[irow] = 0.0;
1412  }
1413  for (jcol = 0; jcol < m_neq; jcol++) {
1414  for (irow = 0; irow < m_neq; irow++) {
1415  m_rowScales[irow] += fabs(*jptr);
1416  jptr++;
1417  }
1418  }
1419 
1420  jptr = &(*(jac.begin()));
1421  for (jcol = 0; jcol < m_neq; jcol++) {
1422  for (irow = 0; irow < m_neq; irow++) {
1423  *jptr /= m_rowScales[irow];
1424  jptr++;
1425  }
1426  }
1427  }
1428  for (irow = 0; irow < m_neq; irow++) {
1429  delta_y[irow] /= m_rowScales[irow];
1430  }
1431  }
1432 
1433 #ifdef DEBUG_HKM_JAC
1434  bool printJacContributions = false;
1435  if (m_time_step_num > 304) {
1436  printJacContributions = false;
1437  }
1438  int focusRow = 10;
1439  int numRows = 2;
1440  double RRow[2];
1441  bool freshJac = true;
1442  RRow[0] = delta_y[focusRow];
1443  RRow[1] = delta_y[focusRow+1];
1444  double Pcutoff = 1.0E-70;
1445  if (!jac.factored()) {
1446  jacBack = jac;
1447  } else {
1448  freshJac = false;
1449  }
1450 #endif
1451  /*
1452  * Solve the system -> This also involves inverting the
1453  * matrix
1454  */
1455  (void) jac.solve(delta_y);
1456 
1457 
1458  /*
1459  * reverse the column scaling if there was any.
1460  */
1461  if (m_colScaling) {
1462  for (irow = 0; irow < m_neq; irow++) {
1463  delta_y[irow] *= m_colScales[irow];
1464  }
1465  }
1466 
1467 #ifdef DEBUG_HKM_JAC
1468  if (printJacContributions) {
1469  for (int iNum = 0; iNum < numRows; iNum++) {
1470  if (iNum > 0) {
1471  focusRow++;
1472  }
1473  double dsum = 0.0;
1474  vector_fp& Jdata = jacBack.data();
1475  double dRow = Jdata[m_neq * focusRow + focusRow];
1476  printf("\n Details on delta_Y for row %d \n", focusRow);
1477  printf(" Value before = %15.5e, delta = %15.5e,"
1478  "value after = %15.5e\n", y_curr[focusRow],
1479  delta_y[focusRow],
1480  y_curr[focusRow] + delta_y[focusRow]);
1481  if (!freshJac) {
1482  printf(" Old Jacobian\n");
1483  }
1484  printf(" col delta_y aij "
1485  "contrib \n");
1486  printf("--------------------------------------------------"
1487  "---------------------------------------------\n");
1488  printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
1489  focusRow, delta_y[focusRow],
1490  dRow, RRow[iNum] / dRow, RRow[iNum]);
1491  dsum += RRow[iNum] / dRow;
1492  for (int ii = 0; ii < m_neq; ii++) {
1493  if (ii != focusRow) {
1494  double aij = Jdata[m_neq * ii + focusRow];
1495  double contrib = aij * delta_y[ii] * (-1.0) / dRow;
1496  dsum += contrib;
1497  if (fabs(contrib) > Pcutoff) {
1498  printf("%6d %15.5e %15.5e %15.5e\n", ii,
1499  delta_y[ii] , aij, contrib);
1500  }
1501  }
1502  }
1503  printf("--------------------------------------------------"
1504  "---------------------------------------------\n");
1505  printf(" %15.5e %15.5e\n",
1506  delta_y[focusRow], dsum);
1507  }
1508  }
1509 
1510 #endif
1511 
1513 }
1514 
1515 double BEulerInt::boundStep(const double* const y,
1516  const double* const step0, int loglevel)
1517 {
1518  int i, i_lower = -1, ifbd = 0, i_fbd = 0;
1519  double fbound = 1.0, f_lowbounds = 1.0, f_delta_bounds = 1.0;
1520  double ff, y_new, ff_alt;
1521  for (i = 0; i < m_neq; i++) {
1522  y_new = y[i] + step0[i];
1523  if ((y_new < (-0.01 * m_ewt[i])) && y[i] >= 0.0) {
1524  ff = 0.9 * (y[i] / (y[i] - y_new));
1525  if (ff < f_lowbounds) {
1526  f_lowbounds = ff;
1527  i_lower = i;
1528  }
1529  }
1530  /*
1531  * Now do a delta bounds
1532  * Increase variables by a factor of 2 only
1533  * decrease variables by a factor of 5 only
1534  */
1535  ff = 1.0;
1536  if ((fabs(y_new) > 2.0 * fabs(y[i])) &&
1537  (fabs(y_new-y[i]) > m_ewt[i])) {
1538  ff = fabs(y[i]/(y_new - y[i]));
1539  ff_alt = fabs(m_ewt[i] / (y_new - y[i]));
1540  ff = std::max(ff, ff_alt);
1541  ifbd = 1;
1542  }
1543  if ((fabs(5.0 * y_new) < fabs(y[i])) &&
1544  (fabs(y_new - y[i]) > m_ewt[i])) {
1545  ff = y[i]/(y_new-y[i]) * (1.0 - 5.0)/5.0;
1546  ff_alt = fabs(m_ewt[i] / (y_new - y[i]));
1547  ff = std::max(ff, ff_alt);
1548  ifbd = 0;
1549  }
1550  if (ff < f_delta_bounds) {
1551  f_delta_bounds = ff;
1552  i_fbd = ifbd;
1553  }
1554  f_delta_bounds = std::min(f_delta_bounds, ff);
1555  }
1556  fbound = std::min(f_lowbounds, f_delta_bounds);
1557  /*
1558  * Report on any corrections
1559  */
1560  if (loglevel > 1) {
1561  if (fbound != 1.0) {
1562  if (f_lowbounds < f_delta_bounds) {
1563  printf("\t\tboundStep: Variable %d causing lower bounds "
1564  "damping of %g\n",
1565  i_lower, f_lowbounds);
1566  } else {
1567  if (ifbd) {
1568  printf("\t\tboundStep: Decrease of Variable %d causing "
1569  "delta damping of %g\n",
1570  i_fbd, f_delta_bounds);
1571  } else {
1572  printf("\t\tboundStep: Increase of variable %d causing"
1573  "delta damping of %g\n",
1574  i_fbd, f_delta_bounds);
1575  }
1576  }
1577  }
1578  }
1579  return fbound;
1580 }
1581 
1582 int BEulerInt::dampStep(double time_curr, const double* y0,
1583  const double* ydot0, const double* step0,
1584  double* y1, double* ydot1, double* step1,
1585  double& s1, GeneralMatrix& jac,
1586  int& loglevel, bool writetitle,
1587  int& num_backtracks)
1588 {
1589  // Compute the weighted norm of the undamped step size step0
1590  double s0 = soln_error_norm(step0);
1591 
1592  // Compute the multiplier to keep all components in bounds
1593  // A value of one indicates that there is no limitation
1594  // on the current step size in the nonlinear method due to
1595  // bounds constraints (either negative values of delta
1596  // bounds constraints.
1597  double fbound = boundStep(y0, step0, loglevel);
1598 
1599  // if fbound is very small, then y0 is already close to the
1600  // boundary and step0 points out of the allowed domain. In
1601  // this case, the Newton algorithm fails, so return an error
1602  // condition.
1603  if (fbound < 1.e-10) {
1604  if (loglevel > 1) {
1605  printf("\t\t\tdampStep: At limits.\n");
1606  }
1607  return -3;
1608  }
1609 
1610  //--------------------------------------------
1611  // Attempt damped step
1612  //--------------------------------------------
1613 
1614  // damping coefficient starts at 1.0
1615  double damp = 1.0;
1616  int j, m;
1617  double ff;
1618  num_backtracks = 0;
1619  for (m = 0; m < NDAMP; m++) {
1620 
1621  ff = fbound*damp;
1622 
1623  // step the solution by the damped step size
1624  /*
1625  * Whenever we update the solution, we must also always
1626  * update the time derivative.
1627  */
1628  for (j = 0; j < m_neq; j++) {
1629  y1[j] = y0[j] + ff*step0[j];
1630  }
1631  calc_ydot(m_order, y1, ydot1);
1632 
1633  // compute the next undamped step, step1[], that would result
1634  // if y1[] were accepted.
1635 
1636  doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel);
1637 
1638 #ifdef DEBUG_HKM
1639  for (j = 0; j < m_neq; j++) {
1640  checkFinite(step1[j]);
1641  checkFinite(y1[j]);
1642  }
1643 #endif
1644  // compute the weighted norm of step1
1645  s1 = soln_error_norm(step1);
1646 
1647  // write log information
1648  if (loglevel > 3) {
1649  print_solnDelta_norm_contrib((const double*) step0,
1650  "DeltaSolnTrial",
1651  (const double*) step1,
1652  "DeltaSolnTrialTest",
1653  "dampNewt: Important Entries for "
1654  "Weighted Soln Updates:",
1655  y0, y1, ff, 5);
1656  }
1657  if (loglevel > 1) {
1658  printf("\t\t\tdampNewt: s0 = %g, s1 = %g, fbound = %g,"
1659  "damp = %g\n", s0, s1, fbound, damp);
1660  }
1661 #ifdef DEBUG_HKM
1662  if (loglevel > 2) {
1663  if (s1 > 1.00000001 * s0 && s1 > 1.0E-5) {
1664  printf("WARNING: Possible Jacobian Problem "
1665  "-> turning on more debugging for this step!!!\n");
1666  print_solnDelta_norm_contrib((const double*) step0,
1667  "DeltaSolnTrial",
1668  (const double*) step1,
1669  "DeltaSolnTrialTest",
1670  "dampNewt: Important Entries for "
1671  "Weighted Soln Updates:",
1672  y0, y1, ff, 5);
1673  loglevel = 4;
1674  }
1675  }
1676 #endif
1677 
1678  // if the norm of s1 is less than the norm of s0, then
1679  // accept this damping coefficient. Also accept it if this
1680  // step would result in a converged solution. Otherwise,
1681  // decrease the damping coefficient and try again.
1682 
1683  if (s1 < 1.0E-5 || s1 < s0) {
1684  if (loglevel > 2) {
1685  if (s1 > s0) {
1686  if (s1 > 1.0) {
1687  printf("\t\t\tdampStep: current trial step and damping"
1688  " coefficient accepted because test step < 1\n");
1689  printf("\t\t\t s1 = %g, s0 = %g\n", s1, s0);
1690  }
1691  }
1692  }
1693  break;
1694  } else {
1695  if (loglevel > 1) {
1696  printf("\t\t\tdampStep: current step rejected: (s1 = %g > "
1697  "s0 = %g)", s1, s0);
1698  if (m < (NDAMP-1)) {
1699  printf(" Decreasing damping factor and retrying");
1700  } else {
1701  printf(" Giving up!!!");
1702  }
1703  printf("\n");
1704  }
1705  }
1706  num_backtracks++;
1707  damp /= DampFactor;
1708  }
1709 
1710  // If a damping coefficient was found, return 1 if the
1711  // solution after stepping by the damped step would represent
1712  // a converged solution, and return 0 otherwise. If no damping
1713  // coefficient could be found, return -2.
1714  if (m < NDAMP) {
1715  if (s1 > 1.0) {
1716  return 0;
1717  } else {
1718  return 1;
1719  }
1720  } else {
1721  if (s1 < 0.5 && (s0 < 0.5)) {
1722  return 1;
1723  }
1724  if (s1 < 1.0) {
1725  return 0;
1726  }
1727  return -2;
1728  }
1729 }
1730 
1731 int BEulerInt::solve_nonlinear_problem(double* const y_comm,
1732  double* const ydot_comm, double CJ,
1733  double time_curr,
1734  GeneralMatrix& jac,
1735  int& num_newt_its,
1736  int& num_linear_solves,
1737  int& num_backtracks,
1738  int loglevel)
1739 {
1740  int m = 0;
1741  bool forceNewJac = false;
1742  double s1=1.e30;
1743 
1744  vector_fp y_curr(y_comm, y_comm + m_neq);
1745  vector_fp ydot_curr(ydot_comm, ydot_comm + m_neq);
1746  vector_fp stp(m_neq, 0.0);
1747  vector_fp stp1(m_neq, 0.0);
1748  vector_fp y_new(m_neq, 0.0);
1749  vector_fp ydot_new(m_neq, 0.0);
1750 
1751  bool frst = true;
1752  num_newt_its = 0;
1753  num_linear_solves = - m_numTotalLinearSolves;
1754  num_backtracks = 0;
1755  int i_backtracks;
1756 
1757  while (1 > 0) {
1758 
1759  /*
1760  * Increment Newton Solve counter
1761  */
1763  num_newt_its++;
1764 
1765 
1766  if (loglevel > 1) {
1767  printf("\t\tSolve_Nonlinear_Problem: iteration %d:\n",
1768  num_newt_its);
1769  }
1770 
1771  // Check whether the Jacobian should be re-evaluated.
1772 
1773  forceNewJac = true;
1774 
1775  if (forceNewJac) {
1776  if (loglevel > 1) {
1777  printf("\t\t\tGetting a new Jacobian and solving system\n");
1778  }
1779  beuler_jac(jac, &m_resid[0], time_curr, CJ, &y_curr[0], &ydot_curr[0],
1780  num_newt_its);
1781  } else {
1782  if (loglevel > 1) {
1783  printf("\t\t\tSolving system with old Jacobian\n");
1784  }
1785  }
1786 
1787  // compute the undamped Newton step
1788  doNewtonSolve(time_curr, &y_curr[0], &ydot_curr[0], &stp[0], jac, loglevel);
1789 
1790  // damp the Newton step
1791  m = dampStep(time_curr, &y_curr[0], &ydot_curr[0], &stp[0], &y_new[0], &ydot_new[0],
1792  &stp1[0], s1, jac, loglevel, frst, i_backtracks);
1793  frst = false;
1794  num_backtracks += i_backtracks;
1795 
1796  /*
1797  * Impose the minimum number of Newton iterations critera
1798  */
1799  if (num_newt_its < m_min_newt_its) {
1800  if (m == 1) {
1801  m = 0;
1802  }
1803  }
1804  /*
1805  * Impose max Newton iteration
1806  */
1807  if (num_newt_its > 20) {
1808  m = -1;
1809  if (loglevel > 1) {
1810  printf("\t\t\tDampnewton unsuccessful (max newts exceeded) sfinal = %g\n", s1);
1811  }
1812  }
1813 
1814  if (loglevel > 1) {
1815  if (m == 1) {
1816  printf("\t\t\tDampNewton iteration successful, nonlin "
1817  "converged sfinal = %g\n", s1);
1818  } else if (m == 0) {
1819  printf("\t\t\tDampNewton iteration successful, get new"
1820  "direction, sfinal = %g\n", s1);
1821  } else {
1822  printf("\t\t\tDampnewton unsuccessful sfinal = %g\n", s1);
1823  }
1824  }
1825 
1826  bool m_filterIntermediate = false;
1827  if (m_filterIntermediate) {
1828  if (m == 0) {
1829  (void) filterNewStep(time_n, &y_new[0], &ydot_new[0]);
1830  }
1831  }
1832  // Exchange new for curr solutions
1833  if (m == 0 || m == 1) {
1834  y_curr = y_new;
1835  calc_ydot(m_order, &y_curr[0], &ydot_curr[0]);
1836  }
1837 
1838  // convergence
1839  if (m == 1) {
1840  goto done;
1841  }
1842 
1843  // If dampStep fails, first try a new Jacobian if an old
1844  // one was being used. If it was a new Jacobian, then
1845  // return -1 to signify failure.
1846  else if (m < 0) {
1847  goto done;
1848  }
1849  }
1850 
1851 done:
1852  // Copy into the return vectors
1853  copy(y_curr.begin(), y_curr.end(), y_comm);
1854  copy(ydot_curr.begin(), ydot_curr.end(), ydot_comm);
1855  // Increment counters
1856  num_linear_solves += m_numTotalLinearSolves;
1857 
1858  double time_elapsed = 0.0;
1859  if (loglevel > 1) {
1860  if (m == 1) {
1861  printf("\t\tNonlinear problem solved successfully in "
1862  "%d its, time elapsed = %g sec\n",
1863  num_newt_its, time_elapsed);
1864  }
1865  }
1866  return m;
1867 }
1868 
1869 void BEulerInt::print_solnDelta_norm_contrib(const double* const solnDelta0,
1870  const char* const s0,
1871  const double* const solnDelta1,
1872  const char* const s1,
1873  const char* const title,
1874  const double* const y0,
1875  const double* const y1,
1876  double damp,
1877  int num_entries)
1878 {
1879  int i, j, jnum;
1880  bool used;
1881  double dmax0, dmax1, error, rel_norm;
1882  printf("\t\t%s currentDamp = %g\n", title, damp);
1883  printf("\t\t I ysoln %10s ysolnTrial "
1884  "%10s weight relSoln0 relSoln1\n", s0, s1);
1885  vector_int imax(num_entries, -1);
1886  printf("\t\t ");
1887  writeline('-', 90);
1888  for (jnum = 0; jnum < num_entries; jnum++) {
1889  dmax1 = -1.0;
1890  for (i = 0; i < m_neq; i++) {
1891  used = false;
1892  for (j = 0; j < jnum; j++) {
1893  if (imax[j] == i) {
1894  used = true;
1895  }
1896  }
1897  if (!used) {
1898  error = solnDelta0[i] / m_ewt[i];
1899  rel_norm = sqrt(error * error);
1900  error = solnDelta1[i] / m_ewt[i];
1901  rel_norm += sqrt(error * error);
1902  if (rel_norm > dmax1) {
1903  imax[jnum] = i;
1904  dmax1 = rel_norm;
1905  }
1906  }
1907  }
1908  if (imax[jnum] >= 0) {
1909  i = imax[jnum];
1910  error = solnDelta0[i] / m_ewt[i];
1911  dmax0 = sqrt(error * error);
1912  error = solnDelta1[i] / m_ewt[i];
1913  dmax1 = sqrt(error * error);
1914  printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e "
1915  "%12.4e %12.4e %12.4e\n",
1916  i, y0[i], solnDelta0[i], y1[i],
1917  solnDelta1[i], m_ewt[i], dmax0, dmax1);
1918  }
1919  }
1920  printf("\t\t ");
1921  writeline('-', 90);
1922 }
1923 
1924 } // End of namespace Cantera
virtual void setPrintSolnOptions(int printSolnStepInterval, int printSolnNumberToTout, int printSolnFirstSteps=0, bool dumpJacobians=false)
This routine controls when the solution is printed.
Definition: BEulerInt.cpp:124
double m_time_final
Final time.
Definition: BEulerInt.h:515
virtual doublereal * ptrColumn(size_t j)=0
Return a pointer to the top of column j, columns are assumed to be contiguous in memory.
virtual void user_out2(const int ifunc, const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
This function may be used to create output at various points in the execution of an application...
int m_numTotalTruncFails
Total Number of time truncation error failures.
Definition: BEulerInt.h:571
double time_error_norm()
Definition: BEulerInt.cpp:639
Dense, Square (not sparse) matrices.
virtual void setProblemType(int probtype)
Set the problem type.
Definition: BEulerInt.cpp:99
int m_printSolnNumberToTout
Number of evenly spaced printouts of the solution If zero, there is no printout from this option defa...
Definition: BEulerInt.h:488
double getPrintTime(double time_current)
Get the next time to print out.
Definition: BEulerInt.cpp:210
virtual void setIterator(IterType t)
Set the linear iterator.
Definition: BEulerInt.cpp:135
int m_printSolnStepInterval
Step Interval at which to print out the solution default = 1; If set to zero, there is no printout...
Definition: BEulerInt.h:481
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
Definition: checkFinite.cpp:43
virtual int nEquations() const
Return the number of equations in the equation system.
int m_numTotalNewtIts
Number of total Newton iterations.
Definition: BEulerInt.h:562
int m_neq
Number of equations in the ode integrator.
Definition: BEulerInt.h:501
virtual double soln_error_norm(const double *const, bool printLargest=false)
Calculate the solution error norm.
Definition: BEulerInt.cpp:1295
void doNewtonSolve(double, double *, double *, double *, GeneralMatrix &, int)
Compute the undamped Newton step.
Definition: BEulerInt.cpp:1349
int m_failure_counter
Failure Counter -> keeps track of the number of consecutive failures.
Definition: BEulerInt.h:469
BEulerMethodType m_method
MethodType is used to specify how the time step is to be chosen.
Definition: BEulerInt.h:401
virtual void clearFactorFlag()
clear the factored flag
bool m_rowScaling
m_rowScaling is a boolean.
Definition: BEulerInt.h:411
virtual int evalResidNJ(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot, doublereal *const resid, const ResidEval_Type_Enum evalType=Base_ResidEval, const int id_x=-1, const doublereal delta_x=0.0)
Evaluate the residual function.
virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal *const y, const doublereal *const ydot, GeneralMatrix &J, doublereal *const resid)
Calculate an analytical Jacobian and the residual at the current time and values. ...
virtual vector_fp::iterator begin()=0
Return an iterator pointing to the first element.
void beuler_jac(GeneralMatrix &J, double *const f, double, double, double *const, double *const, int)
Definition: BEulerInt.cpp:483
double m_reltol
Relative time truncation error tolerances.
Definition: BEulerInt.h:431
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
double m_t0
Initial time at the start of the integration.
Definition: BEulerInt.h:512
const double DampFactor
Dampfactor is the factor by which the damping factor is reduced by when a reduction in step length is...
Definition: BEulerInt.cpp:1288
virtual void writeSolution(int ievent, const double time, const double deltaT, const int time_step_num, const double *y, const double *ydot)
Write out to a file or to standard output the current solution.
Definition: ResidEval.h:135
Generic matrix.
Definition: GeneralMatrix.h:24
Base residual calculation for the Jacobian calculation.
Definition: ResidJacEval.h:28
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
bool m_matrixConditioning
m_matrixConditioning is a boolean.
Definition: BEulerInt.h:423
int m_min_newt_its
Minimum Number of Newton Iterations per nonlinear step. default = 0.
Definition: BEulerInt.h:472
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
int m_print_flag
Determines the level of printing for each time step.
Definition: BEulerInt.h:547
double subtractRD(double a, double b)
This routine subtracts two numbers for one another.
Definition: BEulerInt.cpp:468
vector_fp m_abstol
Vector of absolute time truncation error tolerance when not uniform for all variables.
Definition: BEulerInt.h:442
IterType m_iter
IterType is used to specify how the nonlinear equations are to be relaxed at each time step...
Definition: BEulerInt.h:394
Wrappers for the function evaluators for Nonlinear solvers and Time steppers.
Definition: ResidJacEval.h:51
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
double delta_t_max
Maximum permissible time step.
Definition: BEulerInt.h:526
virtual int calcDeltaSolnVariables(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const delta_y, const doublereal *const solnWeights=0)
Return a vector of delta y's for calculation of the numerical Jacobian.
int solve_nonlinear_problem(double *const y_comm, double *const ydot_comm, double CJ, double time_curr, GeneralMatrix &jac, int &num_newt_its, int &num_linear_solves, int &num_backtracks, int loglevel)
Solve a nonlinear system.
Definition: BEulerInt.cpp:1731
int m_numTotalLinearSolves
Total number of linear iterations.
Definition: BEulerInt.h:565
GeneralMatrix * tdjac_ptr
Pointer to the Jacobian representing the time dependent problem.
Definition: BEulerInt.h:536
virtual void setSolnWeights()
Set the solution weights.
Definition: BEulerInt.cpp:250
virtual int matrixConditioning(doublereal *const matrix, const int nrows, doublereal *const rhs)
Multiply the matrix by another matrix that leads to better conditioning.
int m_printSolnFirstSteps
Number of initial steps that the solution is printed out. default = 0.
Definition: BEulerInt.h:491
virtual void setColumnScales()
Set the column scaling vector at the current time.
Definition: BEulerInt.cpp:266
int m_numInitialConstantDeltaTSteps
Number of initial time steps to take where the time truncation error tolerances are not checked...
Definition: BEulerInt.h:466
void calc_y_pred(int)
Definition: BEulerInt.cpp:590
void error(const std::string &msg)
Write an error message and quit.
Definition: global.cpp:72
int m_order
Current integration order.
Definition: BEulerInt.h:454
int m_nfe
Number of function evaluations.
Definition: BEulerInt.h:553
A class for full (non-sparse) matrices with Fortran-compatible data storage.
Definition: SquareMatrix.h:26
virtual void initializeRJE(double t0, ResidJacEval &func)
Find the initial conditions for y and ydot.
Definition: BEulerInt.cpp:165
virtual void setTolerances(double reltol, size_t n, double *abstol)
Set or reset the number of equations.
Definition: BEulerInt.cpp:78
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
const int NDAMP
Number of damping steps that are carried out before the solution is deemed a failure.
Definition: BEulerInt.cpp:1289
double m_abstols
Absolute time truncation error tolerances, when uniform for all variables.
Definition: BEulerInt.h:437
double time_step_control(int m_order, double time_error_factor)
Definition: BEulerInt.cpp:693
virtual int evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
Evaluate the time tracking equations, if any.
int dampStep(double, const double *, const double *, const double *, double *, double *, double *, double &, GeneralMatrix &, int &, bool, int &)
Definition: BEulerInt.cpp:1582
int m_max_time_step_attempts
Max time steps allowed.
Definition: BEulerInt.h:461
virtual void calcSolnScales(const doublereal t, const doublereal *const y, const doublereal *const y_old, doublereal *const yScales)
Returns a vector of column scale factors that can be used to column scale Jacobians.
double filterNewStep(double, double *, double *)
Filter a new step.
Definition: BEulerInt.cpp:295
int m_jacFormMethod
m_jacFormMethod determines how a matrix is formed.
Definition: BEulerInt.h:405
int m_numTotalConvFails
Total number of convergence failures.
Definition: BEulerInt.h:568
virtual bool factored() const
true if the current factorization is up to date with the matrix
virtual int nEvals() const
Return the total number of function evaluations.
Definition: BEulerInt.cpp:225
double boundStep(const double *const y, const double *const step0, int loglevel)
Bound the Newton step while relaxing the solution.
Definition: BEulerInt.cpp:1515
bool m_colScaling
m_colScaling is a boolean.
Definition: BEulerInt.h:416
virtual int solve(doublereal *b, size_t nrhs=1, size_t ldb=0)=0
Solves the Ax = b system returning x in the b spot.
IterType
Specifies the method used for iteration.
Definition: Integrator.h:40
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
void internalMalloc()
Internal routine that sets up the fixed length storage based on the size of the problem to solve...
Definition: BEulerInt.cpp:230
bool m_dumpJacobians
Dump Jacobians to disk. default false.
Definition: BEulerInt.h:494
virtual doublereal step(double tout)
Integrate the system of equations.
Definition: BEulerInt.cpp:888
Delta residual calculation for the Jacbobian calculation.
Definition: ResidJacEval.h:30
Newton Iteration.
Definition: Integrator.h:42
int m_nJacEval
Number of Jacobian Evaluations and factorization steps (they are the same)
Definition: BEulerInt.h:559
virtual doublereal filterSolnPrediction(const doublereal t, doublereal *const y)
Filter the solution predictions.
void setNonLinOptions(int min_newt_its=0, bool matrixConditioning=false, bool colScaling=false, bool rowScaling=true)
Set the options for the nonlinear method.
Definition: BEulerInt.cpp:140
void calc_ydot(int, double *, double *)
Definition: BEulerInt.cpp:618
double m_hmax
Maximum step size.
Definition: BEulerInt.h:448
int m_itol
If m_itol =1 then each component has an individual value of atol.
Definition: BEulerInt.h:428
virtual bool evalStoppingCritera(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
Evaluate any stopping criteria other than a final time limit.
Base residual calculation for the time-stepping function.
Definition: ResidJacEval.h:26
vector_fp m_ewt
Error Weights. This is a surprisingly important quantity.
Definition: BEulerInt.h:445
int m_time_step_num
Time step number.
Definition: BEulerInt.h:457
void computeResidWts(GeneralMatrix &jac)
Compute Residual Weights.
Definition: BEulerInt.cpp:271