Cantera  2.0
BEulerInt.cpp
Go to the documentation of this file.
1 /**
2  * @file BEulerInt.cpp
3  *
4  */
5 
6 /*
7  * Copyright 2004 Sandia Corporation. Under the terms of Contract
8  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
9  * retains certain rights in this software.
10  * See file License.txt for licensing information.
11  */
12 
14 
15 #include "cantera/base/mdp_allo.h"
16 #include <iostream>
17 
18 using namespace std;
19 using namespace mdp;
20 
21 #define SAFE_DELETE(a) if (a) { delete (a); a = 0; }
22 
23 
24 /*
25  * Blas routines
26  */
27 extern "C" {
28  extern void dcopy_(int*, double*, int*, double*, int*);
29 }
30 namespace Cantera
31 {
32 
33 //================================================================================================
34 /*
35  * Exception thrown when a BEuler error is encountered. We just call the
36  * Cantera Error handler in the initialization list
37  */
38 BEulerErr::BEulerErr(std::string msg) :
39  CanteraError("BEulerInt", msg)
40 {
41 }
42 
43 //================================================================================================
44 /*
45  * Constructor. Default settings: dense jacobian, no user-supplied
46  * Jacobian function, Newton iteration.
47  */
48 BEulerInt::BEulerInt() :
49  m_iter(Newton_Iter),
50  m_method(BEulerVarStep),
51  m_jacFormMethod(BEULER_JAC_NUM),
52  m_rowScaling(true),
53  m_colScaling(false),
54  m_matrixConditioning(false),
55  m_itol(0),
56  m_reltol(1.e-4),
57  m_abstols(1.e-10),
58  m_abstol(0),
59  m_ewt(0),
60  m_hmax(0.0),
61  m_maxord(0),
62  m_time_step_num(0),
63  m_time_step_attempts(0),
64  m_max_time_step_attempts(11000000),
65  m_numInitialConstantDeltaTSteps(0),
66  m_failure_counter(0),
67  m_min_newt_its(0),
68  m_printSolnStepInterval(1),
69  m_printSolnNumberToTout(1),
70  m_printSolnFirstSteps(0),
71  m_dumpJacobians(false),
72  m_neq(0),
73  m_y_n(0),
74  m_y_nm1(0),
75  m_y_pred_n(0),
76  m_ydot_n(0),
77  m_ydot_nm1(0),
78  m_t0(0.0),
79  m_time_final(0.0),
80  time_n(0.0),
81  time_nm1(0.0),
82  time_nm2(0.0),
83  delta_t_n(0.0),
84  delta_t_nm1(0.0),
85  delta_t_nm2(0.0),
86  delta_t_np1(1.0E-8),
87  delta_t_max(1.0E300),
88  m_resid(0),
89  m_residWts(0),
90  m_wksp(0),
91  m_func(0),
92  m_rowScales(0),
93  m_colScales(0),
94  tdjac_ptr(0),
95  m_print_flag(3),
96  m_nfe(0),
97  m_nJacEval(0),
98  m_numTotalNewtIts(0),
99  m_numTotalLinearSolves(0),
100  m_numTotalConvFails(0),
101  m_numTotalTruncFails(0),
102  num_failures(0)
103 {
104 
105 }
106 //================================================================================================
107 /*
108  * Destructor
109  */
111 {
112  mdp::mdp_safe_free((void**) &m_y_n);
113  mdp::mdp_safe_free((void**) &m_y_nm1);
114  mdp::mdp_safe_free((void**) &m_y_pred_n);
115  mdp::mdp_safe_free((void**) &m_ydot_n);
116  mdp::mdp_safe_free((void**) &m_ydot_nm1);
117  mdp::mdp_safe_free((void**) &m_resid);
118  mdp::mdp_safe_free((void**) &m_residWts);
119  mdp::mdp_safe_free((void**) &m_wksp);
120  mdp::mdp_safe_free((void**) &m_ewt);
121  mdp::mdp_safe_free((void**) &m_abstol);
122  mdp::mdp_safe_free((void**) &m_rowScales);
123  mdp::mdp_safe_free((void**) &m_colScales);
124  SAFE_DELETE(tdjac_ptr);
125 }
126 //================================================================================================
127 void BEulerInt::setTolerances(double reltol, size_t n, double* abstol)
128 {
129  m_itol = 1;
130  if (!m_abstol) {
131  m_abstol = mdp_alloc_dbl_1(m_neq, MDP_DBL_NOINIT);
132  }
133  if (static_cast<int>(n) != m_neq) {
134  printf("ERROR n is wrong\n");
135  exit(-1);
136  }
137  for (int i = 0; i < m_neq; i++) {
138  m_abstol[i] = abstol[i];
139  }
140  m_reltol = reltol;
141 }
142 //================================================================================================
143 void BEulerInt::setTolerances(double reltol, double abstol)
144 {
145  m_itol = 0;
146  m_reltol = reltol;
147  m_abstols = abstol;
148 }
149 //================================================================================================
150 void BEulerInt::setProblemType(int jacFormMethod)
151 {
152  m_jacFormMethod = jacFormMethod;
153 }
154 //================================================================================================
155 void BEulerInt::setMethodBEMT(BEulerMethodType t)
156 {
157  m_method = t;
158 }
159 //================================================================================================
160 void BEulerInt::setMaxStep(doublereal hmax)
161 {
162  m_hmax = hmax;
163 }
164 //================================================================================================
165 void BEulerInt::setMaxNumTimeSteps(int maxNumTimeSteps)
166 {
167  m_max_time_step_attempts = maxNumTimeSteps;
168 }
169 //================================================================================================
170 void BEulerInt::setNumInitialConstantDeltaTSteps(int num)
171 {
173 }
174 //================================================================================================
175 /*
176  *
177  * setPrintSolnOptins():
178  *
179  * This routine controls when the solution is printed
180  *
181  * @param printStepInterval If greater than 0, then the
182  * soln is printed every printStepInterval
183  * steps.
184  *
185  * @param printNumberToTout The solution is printed at
186  * regular invervals a total of
187  * "printNumberToTout" times.
188  *
189  * @param printSolnFirstSteps The solution is printed out
190  * the first "printSolnFirstSteps"
191  * steps. After these steps the other
192  * parameters determine the printing.
193  * default = 0
194  *
195  * @param dumpJacobians Dump jacobians to disk.
196  *
197  * default = false
198  *
199  */
200 void BEulerInt::setPrintSolnOptions(int printSolnStepInterval,
201  int printSolnNumberToTout,
202  int printSolnFirstSteps,
203  bool dumpJacobians)
204 {
205  m_printSolnStepInterval = printSolnStepInterval;
206  m_printSolnNumberToTout = printSolnNumberToTout;
207  m_printSolnFirstSteps = printSolnFirstSteps;
208  m_dumpJacobians = dumpJacobians;
209 }
210 //================================================================================================
212 {
213  m_iter = t;
214 }
215 //================================================================================================
216 /*
217  *
218  * setNonLinOptions()
219  *
220  * Set the options for the nonlinear method
221  *
222  * Defaults are set in the .h file. These are the defaults:
223  * min_newt_its = 0
224  * matrixConditioning = false
225  * colScaling = false
226  * rowScaling = true
227  */
228 void BEulerInt::setNonLinOptions(int min_newt_its, bool matrixConditioning,
229  bool colScaling, bool rowScaling)
230 {
231  m_min_newt_its = min_newt_its;
232  m_matrixConditioning = matrixConditioning;
233  m_colScaling = colScaling;
234  m_rowScaling = rowScaling;
235  if (m_colScaling) {
236  if (!m_colScales) {
237  m_colScales = mdp_alloc_dbl_1(m_neq, 1.0);
238  }
239  }
240  if (m_rowScaling) {
241  if (!m_rowScales) {
242  m_rowScales = mdp_alloc_dbl_1(m_neq, 1.0);
243  }
244  }
245 }
246 //================================================================================================
247 /*
248  *
249  * setInitialTimeStep():
250  *
251  * Set the initial time step. Right now, we set the
252  * time step by setting delta_t_np1.
253  */
254 void BEulerInt::setInitialTimeStep(double deltaT)
255 {
256  delta_t_np1 = deltaT;
257 }
258 //================================================================================================
259 /*
260  * setPrintFlag():
261  *
262  */
263 void BEulerInt::setPrintFlag(int print_flag)
264 {
265  m_print_flag = print_flag;
266 }
267 //================================================================================================
268 /*
269  *
270  * initialize():
271  *
272  * Find the initial conditions for y and ydot.
273  */
274 void BEulerInt::initializeRJE(double t0, ResidJacEval& func)
275 {
276  m_neq = func.nEquations();
277  m_t0 = t0;
278  internalMalloc();
279 
280  /*
281  * Get the initial conditions.
282  */
283  func.getInitialConditions(m_t0, m_y_n, m_ydot_n);
284 
285  // Store a pointer to the residual routine in the object
286  m_func = &func;
287 
288  /*
289  * Initialize the various time counters in the object
290  */
291  time_n = t0;
292  time_nm1 = time_n;
293  time_nm2 = time_nm1;
294  delta_t_n = 0.0;
295  delta_t_nm1 = 0.0;
296 }
297 //================================================================================================
298 /*
299  *
300  * reinitialize():
301  *
302  */
303 void BEulerInt::reinitializeRJE(double t0, ResidJacEval& func)
304 {
305  m_neq = func.nEquations();
306  m_t0 = t0;
307  internalMalloc();
308  /*
309  * At the initial time, get the initial conditions and time and store
310  * them into internal storage in the object, my[].
311  */
312  m_t0 = t0;
313  func.getInitialConditions(m_t0, m_y_n, m_ydot_n);
314  /**
315  * Set up the internal weights that are used for testing convergence
316  */
317  setSolnWeights();
318 
319  // Store a pointer to the function
320  m_func = &func;
321 
322 }
323 //================================================================================================
324 /*
325  *
326  * getPrintTime():
327  *
328  */
329 double BEulerInt::getPrintTime(double time_current)
330 {
331  double tnext;
332  if (m_printSolnNumberToTout > 0) {
333  double dt = (m_time_final - m_t0) / m_printSolnNumberToTout;
334  for (int i = 0; i <= m_printSolnNumberToTout; i++) {
335  tnext = m_t0 + dt * i;
336  if (tnext >= time_current) {
337  return tnext;
338  }
339  }
340  }
341  return 1.0E300;
342 }
343 //================================================================================================
344 /*
345  * nEvals():
346  *
347  * Return the total number of function evaluations
348  */
349 int BEulerInt::nEvals() const
350 {
351  return m_nfe;
352 }
353 //================================================================================================
354 /*
355  *
356  * internalMalloc():
357  *
358  * Internal routine that sets up the fixed length storage based on
359  * the size of the problem to solve.
360  */
362 {
363  mdp_realloc_dbl_1(&m_ewt, m_neq, 0, 0.0);
364  mdp_realloc_dbl_1(&m_y_n, m_neq, 0, 0.0);
365  mdp_realloc_dbl_1(&m_y_nm1, m_neq, 0, 0.0);
366  mdp_realloc_dbl_1(&m_y_pred_n, m_neq, 0, 0.0);
367  mdp_realloc_dbl_1(&m_ydot_n, m_neq, 0, 0.0);
368  mdp_realloc_dbl_1(&m_ydot_nm1, m_neq, 0, 0.0);
369  mdp_realloc_dbl_1(&m_resid, m_neq, 0, 0.0);
370  mdp_realloc_dbl_1(&m_residWts, m_neq, 0, 0.0);
371  mdp_realloc_dbl_1(&m_wksp, m_neq, 0, 0.0);
372  if (m_rowScaling) {
373  mdp_realloc_dbl_1(&m_rowScales, m_neq, 0, 1.0);
374  }
375  if (m_colScaling) {
376  mdp_realloc_dbl_1(&m_colScales, m_neq, 0, 1.0);
377  }
379 }
380 //================================================================================================
381 /*
382  * setSolnWeights():
383  *
384  * Set the solution weights
385  * This is a very important routine as it affects quite a few
386  * operations involving convergence.
387  *
388  */
389 void BEulerInt::setSolnWeights()
390 {
391  int i;
392  if (m_itol == 1) {
393  /*
394  * Adjust the atol vector if we are using vector
395  * atol conditions.
396  */
397  // m_func->adjustAtol(m_abstol);
398 
399  for (i = 0; i < m_neq; i++) {
400  m_ewt[i] = m_abstol[i] + m_reltol * 0.5 *
401  (fabs(m_y_n[i]) + fabs(m_y_pred_n[i]));
402  }
403  } else {
404  for (i = 0; i < m_neq; i++) {
405  m_ewt[i] = m_abstols + m_reltol * 0.5 *
406  (fabs(m_y_n[i]) + fabs(m_y_pred_n[i]));
407  }
408  }
409 }
410 //================================================================================================
411 /*
412  *
413  * setColumnScales():
414  *
415  * Set the column scaling vector at the current time
416  */
417 void BEulerInt::setColumnScales()
418 {
419  m_func->calcSolnScales(time_n, m_y_n, m_y_nm1, m_colScales);
420 }
421 //================================================================================================
422 /*
423  * computeResidWts():
424  *
425  * We compute residual weights here, which we define as the L_0 norm
426  * of the Jacobian Matrix, weighted by the solution weights.
427  * This is the proper way to guage the magnitude of residuals. However,
428  * it does need the evaluation of the jacobian, and the implementation
429  * below is slow, but doesn't take up much memory.
430  *
431  * Here a small weighting indicates that the change in solution is
432  * very sensitive to that equation.
433  */
434 void BEulerInt::computeResidWts(GeneralMatrix& jac)
435 {
436  int i, j;
437  double* data = &(*(jac.begin()));
438  double value;
439  for (i = 0; i < m_neq; i++) {
440  m_residWts[i] = fabs(data[i] * m_ewt[0]);
441  for (j = 1; j < m_neq; j++) {
442  value = fabs(data[j*m_neq + i] * m_ewt[j]);
443  m_residWts[i] = std::max(m_residWts[i], value);
444  }
445  }
446 }
447 //================================================================================================
448 /*
449  * filterNewStep():
450  *
451  * void BEulerInt::
452  *
453  */
454 double BEulerInt::filterNewStep(double timeCurrent, double* y_current, double* ydot_current)
455 {
456  return 0.0;
457 }
458 //==================================================================================================
459 static void print_line(const char* str, int n)
460 {
461  for (int i = 0; i < n; i++) {
462  printf("%s", str);
463  }
464  printf("\n");
465 }
466 //==================================================================================================
467 /*
468  * Print out for relevant time step information
469  */
470 static void print_time_step1(int order, int n_time_step, double time,
471  double delta_t_n, double delta_t_nm1,
472  bool step_failed, int num_failures)
473 {
474  const char* string = 0;
475  if (order == 0) {
476  string = "Backward Euler";
477  } else if (order == 1) {
478  string = "Forward/Backward Euler";
479  } else if (order == 2) {
480  string = "Adams-Bashforth/TR";
481  }
482  printf("\n");
483  print_line("=", 80);
484  printf("\nStart of Time Step: %5d Time_n = %9.5g Time_nm1 = %9.5g\n",
485  n_time_step, time, time - delta_t_n);
486  printf("\tIntegration method = %s\n", string);
487  if (step_failed) {
488  printf("\tPreviously attempted step was a failure\n");
489  }
490  if (delta_t_n > delta_t_nm1) {
491  string = "(Increased from previous iteration)";
492  } else if (delta_t_n < delta_t_nm1) {
493  string = "(Decreased from previous iteration)";
494  } else {
495  string = "(same as previous iteration)";
496  }
497  printf("\tdelta_t_n = %8.5e %s", delta_t_n, string);
498  if (num_failures > 0) {
499  printf("\t(Bad_History Failure Counter = %d)", num_failures);
500  }
501  printf("\n\tdelta_t_nm1 = %8.5e\n", delta_t_nm1);
502 }
503 //================================================================================================
504 /*
505  * Print out for relevant time step information
506  */
507 static void print_time_step2(int time_step_num, int order,
508  double time, double time_error_factor,
509  double delta_t_n, double delta_t_np1)
510 {
511  printf("\tTime Step Number %5d was a success: time = %10g\n", time_step_num,
512  time);
513  printf("\t\tEstimated Error\n");
514  printf("\t\t-------------------- = %8.5e\n", time_error_factor);
515  printf("\t\tTolerated Error\n\n");
516  printf("\t- Recommended next delta_t (not counting history) = %g\n",
517  delta_t_np1);
518  printf("\n");
519  print_line("=", 80);
520  printf("\n");
521 }
522 //================================================================================================
523 /*
524  * Print Out descriptive information on why the current step failed
525  */
526 static void print_time_fail(bool convFailure, int time_step_num,
527  double time, double delta_t_n,
528  double delta_t_np1, double time_error_factor)
529 {
530  printf("\n");
531  print_line("=", 80);
532  if (convFailure) {
533  printf("\tTime Step Number %5d experienced a convergence "
534  "failure\n", time_step_num);
535  printf("\tin the non-linear or linear solver\n");
536  printf("\t\tValue of time at failed step = %g\n", time);
537  printf("\t\tdelta_t of the failed step = %g\n",
538  delta_t_n);
539  printf("\t\tSuggested value of delta_t to try next = %g\n",
540  delta_t_np1);
541  } else {
542  printf("\tTime Step Number %5d experienced a truncation error "
543  "failure!\n", time_step_num);
544  printf("\t\tValue of time at failed step = %g\n", time);
545  printf("\t\tdelta_t of the failed step = %g\n",
546  delta_t_n);
547  printf("\t\tSuggested value of delta_t to try next = %g\n",
548  delta_t_np1);
549  printf("\t\tCalculated truncation error factor = %g\n",
550  time_error_factor);
551  }
552  printf("\n");
553  print_line("=", 80);
554 }
555 //================================================================================================
556 /*
557  * Print out the final results and counters
558  */
559 static void print_final(double time, int step_failed,
560  int time_step_num, int num_newt_its,
561  int total_linear_solves, int numConvFails,
562  int numTruncFails, int nfe, int nJacEval)
563 {
564  printf("\n");
565  print_line("=", 80);
566  printf("TIME INTEGRATION ROUTINE HAS FINISHED: ");
567  if (step_failed) {
568  printf(" IT WAS A FAILURE\n");
569  } else {
570  printf(" IT WAS A SUCCESS\n");
571  }
572  printf("\tEnding time = %g\n", time);
573  printf("\tNumber of time steps = %d\n", time_step_num);
574  printf("\tNumber of newt its = %d\n", num_newt_its);
575  printf("\tNumber of linear solves = %d\n", total_linear_solves);
576  printf("\tNumber of convergence failures= %d\n", numConvFails);
577  printf("\tNumber of TimeTruncErr fails = %d\n", numTruncFails);
578  printf("\tNumber of Function evals = %d\n", nfe);
579  printf("\tNumber of Jacobian evals/solvs= %d\n", nJacEval);
580  printf("\n");
581  print_line("=", 80);
582 }
583 //================================================================================================
584 /*
585  * Header info for one line comment about a time step
586  */
587 static void print_lvl1_Header(int nTimes)
588 {
589  printf("\n");
590  if (nTimes) {
591  print_line("-", 80);
592  }
593  printf("time Time Time Time ");
594  if (nTimes == 0) {
595  printf(" START");
596  } else {
597  printf(" (continued)");
598  }
599  printf("\n");
600 
601  printf("step (sec) step Newt Aztc bktr trunc ");
602  printf("\n");
603 
604  printf(" No. Rslt size Its Its stps error |");
605  printf(" comment");
606  printf("\n");
607  print_line("-", 80);
608 }
609 //================================================================================================
610 /*
611  * One line entry about time step
612  * rslt -> 4 letter code
613  */
614 static void print_lvl1_summary(
615  int time_step_num, double time, const char* rslt, double delta_t_n,
616  int newt_its, int aztec_its, int bktr_stps, double time_error_factor,
617  const char* comment)
618 {
619  printf("%6d %11.6g %4s %10.4g %4d %4d %4d %11.4g",
620  time_step_num, time, rslt, delta_t_n, newt_its, aztec_its,
621  bktr_stps, time_error_factor);
622  if (comment) {
623  printf(" | %s", comment);
624  }
625  printf("\n");
626 }
627 //================================================================================================
628 /*
629  * subtractRD():
630  * This routine subtracts 2 numbers. If the difference is less
631  * than 1.0E-14 times the magnitude of the smallest number,
632  * then diff returns an exact zero.
633  * It also returns an exact zero if the difference is less than
634  * 1.0E-300.
635  *
636  * returns: a - b
637  *
638  * This routine is used in numerical differencing schemes in order
639  * to avoid roundoff errors resulting in creating Jacobian terms.
640  * Note: This is a slow routine. However, jacobian errors may cause
641  * loss of convergence. Therefore, in practice this routine
642  * has proved cost-effective.
643  */
644 double subtractRD(double a, double b)
645 {
646  double diff = a - b;
647  double d = std::min(fabs(a), fabs(b));
648  d *= 1.0E-14;
649  double ad = fabs(diff);
650  if (ad < 1.0E-300) {
651  diff = 0.0;
652  }
653  if (ad < d) {
654  diff = 0.0;
655  }
656  return diff;
657 }
658 //================================================================================================
659 /*
660  *
661  * Function called by BEuler to evaluate the Jacobian matrix and the
662  * current residual at the current time step.
663  * @param N = The size of the equation system
664  * @param J = Jacobian matrix to be filled in
665  * @param f = Right hand side. This routine returns the current
666  * value of the rhs (output), so that it does
667  * not have to be computed again.
668  *
669  */
670 void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f,
671  double time_curr, double CJ,
672  double* const y,
673  double* const ydot,
674  int num_newt_its)
675 {
676  int i, j;
677  double* col_j;
678  double ysave, ydotsave, dy;
679  /**
680  * Clear the factor flag
681  */
682  J.clearFactorFlag();
683 
684 
685  if (m_jacFormMethod & BEULER_JAC_ANAL) {
686  /********************************************************************
687  * Call the function to get a jacobian.
688  */
689  m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
690 #ifdef DEBUG_HKM
691  //double dddd = J(89, 89);
692  //checkFinite(dddd);
693 #endif
694  m_nJacEval++;
695  m_nfe++;
696  } else {
697  /*******************************************************************
698  * Generic algorithm to calculate a numerical Jacobian
699  */
700  /*
701  * Calculate the current value of the rhs given the
702  * current conditions.
703  */
704 
705  m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
706  m_nfe++;
707  m_nJacEval++;
708 
709 
710  /*
711  * Malloc a vector and call the function object to return a set of
712  * deltaY's that are appropriate for calculating the numerical
713  * derivative.
714  */
715  double* dyVector = mdp::mdp_alloc_dbl_1(m_neq, MDP_DBL_NOINIT);
716  m_func->calcDeltaSolnVariables(time_curr, y, m_y_nm1, dyVector,
717  m_ewt);
718 #ifdef DEBUG_HKM
719  bool print_NumJac = false;
720  if (print_NumJac) {
721  FILE* idy = fopen("NumJac.csv", "w");
722  fprintf(idy, "Unk m_ewt y "
723  "dyVector ResN\n");
724  for (int iii = 0; iii < m_neq; iii++) {
725  fprintf(idy, " %4d %16.8e %16.8e %16.8e %16.8e \n",
726  iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
727  }
728  fclose(idy);
729  }
730 #endif
731  /*
732  * Loop over the variables, formulating a numerical derivative
733  * of the dense matrix.
734  * For the delta in the variable, we will use a variety of approaches
735  * The original approach was to use the error tolerance amount.
736  * This may not be the best approach, as it could be overly large in
737  * some instances and overly small in others.
738  * We will first protect from being overly small, by using the usual
739  * sqrt of machine precision approach, i.e., 1.0E-7,
740  * to bound the lower limit of the delta.
741  */
742  for (j = 0; j < m_neq; j++) {
743 
744 
745  /*
746  * Get a pointer into the column of the matrix
747  */
748 
749 
750  col_j = (double*) J.ptrColumn(j);
751  ysave = y[j];
752  dy = dyVector[j];
753  //dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
754 
755  y[j] = ysave + dy;
756  dy = y[j] - ysave;
757  ydotsave = ydot[j];
758  ydot[j] += dy * CJ;
759  /*
760  * Call the functon
761  */
762 
763 
764  m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, m_wksp,
765  JacDelta_ResidEval, j, dy);
766  m_nfe++;
767  double diff;
768  for (i = 0; i < m_neq; i++) {
769  diff = subtractRD(m_wksp[i], f[i]);
770  col_j[i] = diff / dy;
771  //col_j[i] = (m_wksp[i] - f[i])/dy;
772  }
773 
774  y[j] = ysave;
775  ydot[j] = ydotsave;
776 
777  }
778  /*
779  * Release memory
780  */
781  mdp::mdp_safe_free((void**) &dyVector);
782  }
783 
784 
785 }
786 
787 
788 /*
789  * Function to calculate the predicted solution vector, m_y_pred_n for the
790  * (n+1)th time step. This routine can be used by a first order - forward
791  * Euler / backward Euler predictor / corrector method or for a second order
792  * Adams-Bashforth / Trapezoidal Rule predictor / corrector method. See Nachos
793  * documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL - 83282 for
794  * more information.
795  *
796  * variables:
797  *
798  * on input:
799  *
800  * N - number of unknowns
801  * order - indicates order of method
802  * = 1 -> first order forward Euler/backward Euler
803  * predictor/corrector
804  * = 2 -> second order Adams-Bashforth/Trapezoidal Rule
805  * predictor/corrector
806  *
807  * delta_t_n - magnitude of time step at time n (i.e., = t_n+1 - t_n)
808  * delta_t_nm1 - magnitude of time step at time n - 1 (i.e., = t_n - t_n-1)
809  * y_n[] - solution vector at time n
810  * y_dot_n[] - acceleration vector from the predictor at time n
811  * y_dot_nm1[] - acceleration vector from the predictor at time n - 1
812  *
813  * on output:
814  *
815  * m_y_pred_n[] - predicted solution vector at time n + 1
816  */
817 void BEulerInt::calc_y_pred(int order)
818 {
819  int i;
820  double c1, c2;
821  switch (order) {
822  case 0:
823  case 1:
824  c1 = delta_t_n;
825  for (i = 0; i < m_neq; i++) {
826  m_y_pred_n[i] = m_y_n[i] + c1 * m_ydot_n[i];
827  }
828  break;
829  case 2:
830  c1 = delta_t_n * (2.0 + delta_t_n / delta_t_nm1) / 2.0;
831  c2 = (delta_t_n * delta_t_n) / (delta_t_nm1 * 2.0);
832  for (i = 0; i < m_neq; i++) {
833  m_y_pred_n[i] = m_y_n[i] + c1 * m_ydot_n[i] - c2 * m_ydot_nm1[i];
834  }
835  break;
836  }
837 
838  /*
839  * Filter the predictions.
840  */
841  m_func->filterSolnPrediction(time_n, m_y_pred_n);
842 
843 } /* calc_y_pred */
844 
845 
846 /* Function to calculate the acceleration vector ydot for the first or
847  * second order predictor/corrector time integrator. This routine can be
848  * called by a first order - forward Euler / backward Euler predictor /
849  * corrector or for a second order Adams - Bashforth / Trapezoidal Rule
850  * predictor / corrector. See Nachos documentation Sand86-1816 and Gresho,
851  * Lee, Sani LLNL report UCRL - 83282 for more information.
852  *
853  * variables:
854  *
855  * on input:
856  *
857  * N - number of local unknowns on the processor
858  * This is equal to internal plus border unknowns.
859  * order - indicates order of method
860  * = 1 -> first order forward Euler/backward Euler
861  * predictor/corrector
862  * = 2 -> second order Adams-Bashforth/Trapezoidal Rule
863  * predictor/corrector
864  *
865  * delta_t_n - Magnitude of the current time step at time n
866  * (i.e., = t_n - t_n-1)
867  * y_curr[] - Current Solution vector at time n
868  * y_nm1[] - Solution vector at time n-1
869  * ydot_nm1[] - Acceleration vector at time n-1
870  *
871  * on output:
872  *
873  * ydot_curr[] - Current acceleration vector at time n
874  *
875  * Note we use the current attribute to denote the possibility that
876  * y_curr[] may not be equal to m_y_n[] during the nonlinear solve
877  * because we may be using a look-ahead scheme.
878  */
879 void BEulerInt::
880 calc_ydot(int order, double* y_curr, double* ydot_curr)
881 {
882  int i;
883  double c1;
884  switch (order) {
885  case 0:
886  case 1: /* First order forward Euler/backward Euler */
887  c1 = 1.0 / delta_t_n;
888  for (i = 0; i < m_neq; i++) {
889  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]);
890  }
891  return;
892  case 2: /* Second order Adams-Bashforth / Trapezoidal Rule */
893  c1 = 2.0 / delta_t_n;
894  for (i = 0; i < m_neq; i++) {
895  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i];
896  }
897  return;
898  }
899 } /************* END calc_ydot () ****************************************/
900 
901 /* This function calculates the time step truncation error estimate
902  * from a very simple formula based on Gresho et al. This routine can be
903  * called for a
904  * first order - forward Euler/backward Euler predictor/ corrector and
905  * for a
906  * second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
907  * Nachos documentation Sand86-1816 and Gresho, Lee, LLNL report
908  * UCRL - 83282
909  * for more information.
910  *
911  * variables:
912  *
913  * on input:
914  *
915  * abs_error - Generic absolute error tolerance
916  * rel_error - Generic realtive error tolerance
917  * x_coor[] - Solution vector from the implicit corrector
918  * x_pred_n[] - Solution vector from the explicit predictor
919  *
920  * on output:
921  *
922  * delta_t_n - Magnitude of next time step at time t_n+1
923  * delta_t_nm1 - Magnitude of previous time step at time t_n
924  */
926 {
927  int i;
928  double rel_norm, error;
929 #ifdef DEBUG_HKM
930 #define NUM_ENTRIES 5
931  if (m_print_flag > 2) {
932  int imax[NUM_ENTRIES], j, jnum;
933  double dmax;
934  bool used;
935  printf("\t\ttime step truncation error contributors:\n");
936  printf("\t\t I entry actual predicted "
937  " weight ydot\n");
938  printf("\t\t");
939  print_line("-", 70);
940  for (j = 0; j < NUM_ENTRIES; j++) {
941  imax[j] = -1;
942  }
943  for (jnum = 0; jnum < NUM_ENTRIES; jnum++) {
944  dmax = -1.0;
945  for (i = 0; i < m_neq; i++) {
946  used = false;
947  for (j = 0; j < jnum; j++) {
948  if (imax[j] == i) {
949  used = true;
950  }
951  }
952  if (!used) {
953  error = (m_y_n[i] - m_y_pred_n[i]) / m_ewt[i];
954  rel_norm = sqrt(error * error);
955  if (rel_norm > dmax) {
956  imax[jnum] = i;
957  dmax = rel_norm;
958  }
959  }
960  }
961  if (imax[jnum] >= 0) {
962  i = imax[jnum];
963  printf("\t\t%4d %12.4e %12.4e %12.4e %12.4e %12.4e\n",
964  i, dmax, m_y_n[i], m_y_pred_n[i], m_ewt[i], m_ydot_n[i]);
965  }
966  }
967  printf("\t\t");
968  print_line("-", 70);
969  }
970 #endif
971  rel_norm = 0.0;
972  for (i = 0; i < m_neq; i++) {
973  error = (m_y_n[i] - m_y_pred_n[i]) / m_ewt[i];
974  rel_norm += (error * error);
975  }
976  rel_norm = sqrt(rel_norm / m_neq);
977  return rel_norm;
978 }
979 
980 /*************************************************************************
981  * Time step control function for the selection of the time step size based on
982  * a desired accuracy of time integration and on an estimate of the relative
983  * error of the time integration process. This routine can be called for a
984  * first order - forward Euler/backward Euler predictor/ corrector and for a
985  * second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
986  * Nachos documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL -
987  * 83282 for more information.
988  *
989  * variables:
990  *
991  * on input:
992  *
993  * order - indicates order of method
994  * = 1 -> first order forward Euler/backward Euler
995  * predictor/corrector
996  * = 2 -> second order forward Adams-Bashforth/Trapezoidal
997  * rule predictor/corrector
998  *
999  * delta_t_n - Magnitude of time step at time t_n
1000  * delta_t_nm1 - Magnitude of time step at time t_n-1
1001  * rel_error - Generic realtive error tolerance
1002  * time_error_factor - Estimated value of the time step truncation error
1003  * factor. This value is a ratio of the computed
1004  * error norms. The premultiplying constants
1005  * and the power are not yet applied to normalize the
1006  * predictor/corrector ratio. (see output value)
1007  *
1008  * on output:
1009  *
1010  * return - delta_t for the next time step
1011  * If delta_t is negative, then the current time step is
1012  * rejected because the time-step truncation error is
1013  * too large. The return value will contain the negative
1014  * of the recommended next time step.
1015  *
1016  * time_error_factor - This output value is normalized so that
1017  * values greater than one indicate the current time
1018  * integration error is greater than the user
1019  * specified magnitude.
1020  */
1021 double BEulerInt::time_step_control(int order, double time_error_factor)
1022 {
1023  double factor = 0.0, power = 0.0, delta_t;
1024  const char* yo = "time_step_control";
1025 
1026  /*
1027  * Special case time_error_factor so that zeroes don't cause a problem.
1028  */
1029  time_error_factor = std::max(1.0E-50, time_error_factor);
1030 
1031  /*
1032  * Calculate the factor for the change in magnitude of time step.
1033  */
1034  switch (order) {
1035  case 1:
1036  factor = 1.0/(2.0 *(time_error_factor));
1037  power = 0.5;
1038  break;
1039  case 2:
1040  factor = 1.0/(3.0 * (1.0 + delta_t_nm1 / delta_t_n)
1041  * (time_error_factor));
1042  power = 0.3333333333333333;
1043  }
1044  factor = pow(factor, power);
1045  if (factor < 0.5) {
1046  if (m_print_flag > 1) {
1047  printf("\t%s: WARNING - Current time step will be chucked\n", yo);
1048  printf("\t\tdue to a time step truncation error failure.\n");
1049  }
1050  delta_t = - 0.5 * delta_t_n;
1051  } else {
1052  factor = std::min(factor, 1.5);
1053  delta_t = factor * delta_t_n;
1054  }
1055  return delta_t;
1056 } /************ END of time_step_control()********************************/
1057 //================================================================================================
1058 /**************************************************************************
1059  *
1060  * integrate():
1061  *
1062  * defaults are located in the .h file. They are as follows:
1063  * time_init = 0.0
1064  */
1065 double BEulerInt::integrateRJE(double tout, double time_init)
1066 {
1067  double time_current;
1068  bool weAreNotFinished = true;
1069  m_time_final = tout;
1070  int flag = SUCCESS;
1071  /**
1072  * Initialize the time step number to zero. step will increment so that
1073  * the first time step is number 1
1074  */
1075  m_time_step_num = 0;
1076 
1077 
1078  /*
1079  * Do the integration a step at a time
1080  */
1081  int istep = 0;
1082  int printStep = 0;
1083  bool doPrintSoln = false;
1084  time_current = time_init;
1085  time_n = time_init;
1086  time_nm1 = time_init;
1087  time_nm2 = time_init;
1088  m_func->evalTimeTrackingEqns(time_current, 0.0, m_y_n, m_ydot_n);
1089  double print_time = getPrintTime(time_current);
1090  if (print_time == time_current) {
1091  m_func->writeSolution(4, time_current, delta_t_n,
1092  istep, m_y_n, m_ydot_n);
1093  }
1094  /*
1095  * We print out column headers here for the case of
1096  */
1097  if (m_print_flag == 1) {
1098  print_lvl1_Header(0);
1099  }
1100  /*
1101  * Call a different user routine at the end of each step,
1102  * that will probably print to a file.
1103  */
1104  m_func->user_out2(0, time_current, 0.0, m_y_n, m_ydot_n);
1105 
1106  do {
1107 
1108  print_time = getPrintTime(time_current);
1109  if (print_time >= tout) {
1110  print_time = tout;
1111  }
1112 
1113  /************************************************************
1114  * Step the solution
1115  */
1116  time_current = step(tout);
1117  istep++;
1118  printStep++;
1119  /***********************************************************/
1120  if (time_current < 0.0) {
1121  if (time_current == -1234.) {
1122  time_current = 0.0;
1123  } else {
1124  time_current = -time_current;
1125  }
1126  flag = FAILURE;
1127  }
1128 
1129  if (flag != FAILURE) {
1130  bool retn =
1131  m_func->evalStoppingCritera(time_current, delta_t_n,
1132  m_y_n, m_ydot_n);
1133  if (retn) {
1134  weAreNotFinished = false;
1135  doPrintSoln = true;
1136  }
1137  }
1138 
1139  /*
1140  * determine conditional printing of soln
1141  */
1142  if (time_current >= print_time) {
1143  doPrintSoln = true;
1144  }
1145  if (m_printSolnStepInterval == printStep) {
1146  doPrintSoln = true;
1147  }
1148  if (m_printSolnFirstSteps > istep) {
1149  doPrintSoln = true;
1150  }
1151 
1152  /*
1153  * Evaluate time integrated quantities that are calculated at the
1154  * end of every successful time step.
1155  */
1156  if (flag != FAILURE) {
1157  m_func->evalTimeTrackingEqns(time_current, delta_t_n,
1158  m_y_n, m_ydot_n);
1159  }
1160 
1161  /*
1162  * Call the printout routine.
1163  */
1164  if (doPrintSoln) {
1165  m_func->writeSolution(1, time_current, delta_t_n,
1166  istep, m_y_n, m_ydot_n);
1167  printStep = 0;
1168  doPrintSoln = false;
1169  if (m_print_flag == 1) {
1170  print_lvl1_Header(1);
1171  }
1172  }
1173  /*
1174  * Call a different user routine at the end of each step,
1175  * that will probably print to a file.
1176  */
1177  if (flag == FAILURE) {
1178  m_func->user_out2(-1, time_current, delta_t_n, m_y_n, m_ydot_n);
1179  } else {
1180  m_func->user_out2(1, time_current, delta_t_n, m_y_n, m_ydot_n);
1181  }
1182 
1183  } while (time_current < tout &&
1184  m_time_step_attempts < m_max_time_step_attempts &&
1185  flag == SUCCESS && weAreNotFinished);
1186 
1187  /*
1188  * Check current time against the max solution time.
1189  */
1190  if (time_current >= tout) {
1191  printf("Simulation completed time integration in %d time steps\n",
1192  m_time_step_num);
1193  printf("Final Time: %e\n\n", time_current);
1194  } else if (m_time_step_attempts >= m_max_time_step_attempts) {
1195  printf("Simulation ran into time step attempt limit in"
1196  "%d time steps\n",
1197  m_time_step_num);
1198  printf("Final Time: %e\n\n", time_current);
1199  } else if (flag == FAILURE) {
1200  printf("ERROR: time stepper failed at time = %g\n", time_current);
1201  }
1202 
1203  /*
1204  * Print out the final results and counters.
1205  */
1206  print_final(time_n, flag, m_time_step_num, m_numTotalNewtIts,
1209 
1210  /*
1211  * Call a different user routine at the end of each step,
1212  * that will probably print to a file.
1213  */
1214  m_func->user_out2(2, time_current, delta_t_n, m_y_n, m_ydot_n);
1215 
1216 
1217  if (flag != SUCCESS) {
1218  throw BEulerErr(" BEuler error encountered.");
1219  }
1220  return time_current;
1221 }
1222 
1223 /**************************************************************************
1224  *
1225  * step():
1226  *
1227  * This routine advances the calculations one step using a predictor
1228  * corrector approach. We use an implicit algorithm here.
1229  *
1230  */
1231 double BEulerInt::step(double t_max)
1232 {
1233  double CJ;
1234  int one = 1;
1235  bool step_failed = false;
1236  bool giveUp = false;
1237  bool convFailure = false;
1238  const char* rslt;
1239  double time_error_factor = 0.0;
1240  double normFilter = 0.0;
1241  int numTSFailures = 0;
1242  int bktr_stps = 0;
1243  int nonlinearloglevel = m_print_flag;
1244  int num_newt_its = 0;
1245  int aztec_its = 0;
1246  string comment;
1247  /*
1248  * Increment the time counter - May have to be taken back,
1249  * if time step is found to be faulty.
1250  */
1251  m_time_step_num++;
1252 
1253  /**
1254  * Loop here until we achieve a successful step or we set the giveUp
1255  * flag indicating that repeated errors have occurred.
1256  */
1257  do {
1258  m_time_step_attempts++;
1259  comment.clear();
1260 
1261  /*
1262  * Possibly adjust the delta_t_n value for this time step from the
1263  * recommended delta_t_np1 value determined in the previous step
1264  * due to maximum time step constraints or other occurences,
1265  * known to happen at a given time.
1266  */
1267  if ((time_n + delta_t_np1) >= t_max) {
1268  delta_t_np1 =t_max - time_n;
1269  }
1270 
1271  if (delta_t_np1 >= delta_t_max) {
1272  delta_t_np1 = delta_t_max;
1273  }
1274 
1275  /*
1276  * Increment the delta_t counters and the time for the current
1277  * time step.
1278  */
1279 
1280  delta_t_nm2 = delta_t_nm1;
1281  delta_t_nm1 = delta_t_n;
1282  delta_t_n = delta_t_np1;
1283  time_n += delta_t_n;
1284 
1285  /*
1286  * Determine the integration order of the current step.
1287  *
1288  * Special case for start-up of time integration procedure
1289  * First time step = Do a predictor step as we
1290  * have recently added an initial
1291  * ydot input option. And, setting ydot=0
1292  * is equivalent to not doing a
1293  * predictor step.
1294  * Second step = If 2nd order method, do a first order
1295  * step for this time-step, only.
1296  *
1297  * If 2nd order method with a constant time step, the
1298  * first and second steps are 1/10 the specified step, and
1299  * the third step is 8/10 the specified step. This reduces
1300  * the error asociated with using lower order
1301  * integration on the first two steps. (RCS 11-6-97)
1302  *
1303  * If the previous time step failed for one reason or another,
1304  * do a linear step. It's more robust.
1305  */
1306  if (m_time_step_num == 1) {
1307  m_order = 1; /* Backward Euler */
1308  } else if (m_time_step_num == 2) {
1309  m_order = 1; /* Forward/Backward Euler */
1310  } else if (step_failed) {
1311  m_order = 1; /* Forward/Backward Euler */
1312  } else if (m_time_step_num > 2) {
1313  m_order = 1; /* Specified
1314  Predictor/Corrector
1315  - not implemented */
1316  }
1317 
1318  /*
1319  * Print out an initial statement about the step.
1320  */
1321  if (m_print_flag > 1) {
1322  print_time_step1(m_order, m_time_step_num, time_n, delta_t_n,
1323  delta_t_nm1, step_failed, m_failure_counter);
1324  }
1325 
1326  /*
1327  * Calculate the predicted solution, m_y_pred_n, for the current
1328  * time step.
1329  */
1331 
1332  /*
1333  * HKM - Commented this out. I may need it for particles later.
1334  * If Solution bounds checking is turned on, we need to crop the
1335  * predicted solution to make sure bounds are enforced
1336  *
1337  *
1338  * cropNorm = 0.0;
1339  * if (Cur_Realm->Realm_Nonlinear.Constraint_Backtracking_Flag ==
1340  * Constraint_Backtrack_Enable) {
1341  * cropNorm = cropPredictor(mesh, x_pred_n, abs_time_error,
1342  * m_reltol);
1343  */
1344 
1345  /*
1346  * Save the old solution, before overwriting with the new solution
1347  * - use
1348  */
1349  mdp_copy_dbl_1(m_y_nm1, m_y_n, m_neq);
1350 
1351  /*
1352  * Use the predicted value as the initial guess for the corrector
1353  * loop, for
1354  * every step other than the first step.
1355  */
1356  if (m_order > 0) {
1357  mdp_copy_dbl_1(m_y_n, m_y_pred_n, m_neq);
1358  }
1359 
1360  /*
1361  * Save the old time derivative, if necessary, before it is
1362  * overwritten.
1363  * This overwrites ydot_nm1, losing information from the previous time
1364  * step.
1365  */
1366  mdp_copy_dbl_1(m_ydot_nm1, m_ydot_n, m_neq);
1367 
1368  /*
1369  * Calculate the new time derivative, ydot_n, that is consistent
1370  * with the
1371  * initial guess for the corrected solution vector.
1372  *
1373  */
1374  calc_ydot(m_order, m_y_n, m_ydot_n);
1375 
1376  /*
1377  * Calculate CJ, the coefficient for the jacobian corresponding to the
1378  * derivative of the residual wrt to the acceleration vector.
1379  */
1380  if (m_order < 2) {
1381  CJ = 1.0 / delta_t_n;
1382  } else {
1383  CJ = 2.0 / delta_t_n;
1384  }
1385 
1386  /*
1387  * Calculate a new Solution Error Weighting vector
1388  */
1389  setSolnWeights();
1390 
1391  /*
1392  * Solve the system of equations at the current time step.
1393  * Note - x_corr_n and x_dot_n are considered to be updated,
1394  * on return from this solution.
1395  */
1396  int ierror = solve_nonlinear_problem(m_y_n, m_ydot_n,
1397  CJ, time_n, *tdjac_ptr, num_newt_its,
1398  aztec_its, bktr_stps,
1399  nonlinearloglevel);
1400  /*
1401  * Set the appropriate flags if a convergence failure is detected.
1402  */
1403  if (ierror < 0) { /* Step failed */
1404  convFailure = true;
1405  step_failed = true;
1406  rslt = "fail";
1408  m_failure_counter +=3;
1409  if (m_print_flag > 1) {
1410  printf("\tStep is Rejected, nonlinear problem didn't converge,"
1411  "ierror = %d\n", ierror);
1412  }
1413  } else { /* Step succeeded */
1414  convFailure = false;
1415  step_failed = false;
1416  rslt = "done";
1417 
1418  /*
1419  * Apply a filter to a new successful step
1420  */
1421  normFilter = filterNewStep(time_n, m_y_n, m_ydot_n);
1422  if (normFilter > 1.0) {
1423  convFailure = true;
1424  step_failed = true;
1425  rslt = "filt";
1426  if (m_print_flag > 1) {
1427  printf("\tStep is Rejected, too large filter adjustment = %g\n",
1428  normFilter);
1429  }
1430  } else if (normFilter > 0.0) {
1431  if (normFilter > 0.3) {
1432  if (m_print_flag > 1) {
1433  printf("\tStep was filtered, norm = %g, next "
1434  "time step adjusted\n", normFilter);
1435  }
1436  } else {
1437  if (m_print_flag > 1) {
1438  printf("\tStep was filtered, norm = %g\n", normFilter);
1439  }
1440  }
1441  }
1442  }
1443 
1444  /*
1445  * Calculate the time step truncation error for the current step.
1446  */
1447  if (!step_failed) {
1448  time_error_factor = time_error_norm();
1449  } else {
1450  time_error_factor = 1000.;
1451  }
1452 
1453  /*
1454  * Dynamic time step control- delta_t_n, delta_t_nm1 are set here.
1455  */
1456  if (step_failed) {
1457  /*
1458  * For convergence failures, decrease the step-size by a factor of
1459  * 4 and try again.
1460  */
1461  delta_t_np1 = 0.25 * delta_t_n;
1462  } else if (m_method == BEulerVarStep) {
1463 
1464  /*
1465  * If we are doing a predictor/corrector method, and we are
1466  * past a certain number of time steps given by the input file
1467  * then either correct the DeltaT for the next time step or
1468  *
1469  */
1470  if ((m_order > 0) &&
1472  delta_t_np1 = time_step_control(m_order, time_error_factor);
1473  if (normFilter > 0.1) {
1474  if (delta_t_np1 > delta_t_n) {
1475  delta_t_np1 = delta_t_n;
1476  }
1477  }
1478 
1479  /*
1480  * Check for Current time step failing due to violation of
1481  * time step
1482  * truncation bounds.
1483  */
1484  if (delta_t_np1 < 0.0) {
1486  step_failed = true;
1487  delta_t_np1 = -delta_t_np1;
1488  m_failure_counter += 2;
1489  comment += "TIME TRUNC FAILURE";
1490  rslt = "TRNC";
1491  }
1492 
1493  /*
1494  * Prevent churning of the time step by not increasing the
1495  * time step,
1496  * if the recent "History" of the time step behavior is still bad
1497  */
1498  else if (m_failure_counter > 0) {
1499  delta_t_np1 = std::min(delta_t_np1, delta_t_n);
1500  }
1501  } else {
1502  delta_t_np1 = delta_t_n;
1503  }
1504 
1505  /* Decrease time step if a lot of Newton Iterations are
1506  * taken.
1507  * The idea being if more or less Newton iteration are taken
1508  * than the target number of iterations, then adjust the time
1509  * step downwards so that the target number of iterations or lower
1510  * is achieved. This
1511  * should prevent step failure by too many Newton iterations because
1512  * the time step becomes too large. CCO
1513  * hkm -> put in num_new_its min of 3 because the time step
1514  * was being altered even when num_newt_its == 1
1515  */
1516  int max_Newton_steps = 10000;
1517  int target_num_iter = 5;
1518  if (num_newt_its > 3000 && !step_failed) {
1519  if (max_Newton_steps != target_num_iter) {
1520  double iter_diff = num_newt_its - target_num_iter;
1521  double iter_adjust_zone = max_Newton_steps - target_num_iter;
1522  double target_time_step = delta_t_n
1523  *(1.0 - iter_diff*fabs(iter_diff)/
1524  ((2.0*iter_adjust_zone*iter_adjust_zone)));
1525  target_time_step = std::max(0.5*delta_t_n, target_time_step);
1526  if (target_time_step < delta_t_np1) {
1527  printf("\tNext time step will be decreased from %g to %g"
1528  " because of new its restraint\n",
1529  delta_t_np1, target_time_step);
1530  delta_t_np1 = target_time_step;
1531  }
1532  }
1533  }
1534 
1535 
1536  }
1537 
1538  /*
1539  * The final loop in the time stepping algorithm depends on whether the
1540  * current step was a success or not.
1541  */
1542  if (step_failed) {
1543  /*
1544  * Increment the counter indicating the number of consecutive
1545  * failures
1546  */
1547  numTSFailures++;
1548  /*
1549  * Print out a statement about the failure of the time step.
1550  */
1551  if (m_print_flag > 1) {
1552  print_time_fail(convFailure, m_time_step_num, time_n, delta_t_n,
1553  delta_t_np1, time_error_factor);
1554  } else if (m_print_flag == 1) {
1555  print_lvl1_summary(m_time_step_num, time_n, rslt, delta_t_n,
1556  num_newt_its, aztec_its, bktr_stps,
1557  time_error_factor,
1558  comment.c_str());
1559  }
1560 
1561  /*
1562  * Change time step counters back to the previous step before
1563  * the failed
1564  * time step occurred.
1565  */
1566  time_n -= delta_t_n;
1567  delta_t_n = delta_t_nm1;
1568  delta_t_nm1 = delta_t_nm2;
1569 
1570  /*
1571  * Replace old solution vector and time derivative solution vector.
1572  */
1573  dcopy_(&m_neq, m_y_nm1, &one, m_y_n, &one);
1574  dcopy_(&m_neq, m_ydot_nm1, &one, m_ydot_n, &one);
1575  /*
1576  * Decide whether to bail on the whole loop
1577  */
1578  if (numTSFailures > 35) {
1579  giveUp = true;
1580  }
1581  }
1582 
1583  /*
1584  * Do processing for a successful step.
1585  */
1586  else {
1587 
1588  /*
1589  * Decrement the number of consequative failure counter.
1590  */
1592 
1593  /*
1594  * Print out final results of a successfull time step.
1595  */
1596  if (m_print_flag > 1) {
1597  print_time_step2(m_time_step_num, m_order, time_n, time_error_factor,
1598  delta_t_n, delta_t_np1);
1599  } else if (m_print_flag == 1) {
1600  print_lvl1_summary(m_time_step_num, time_n, " ", delta_t_n,
1601  num_newt_its, aztec_its, bktr_stps, time_error_factor,
1602  comment.c_str());
1603  }
1604 
1605  /*
1606  * Output information at the end of every successful time step, if
1607  * requested.
1608  *
1609  * fill in
1610  */
1611 
1612 
1613  }
1614  } while (step_failed && !giveUp);
1615 
1616  /*
1617  * Send back the overall result of the time step.
1618  */
1619  if (step_failed) {
1620  if (time_n == 0.0) {
1621  return -1234.0;
1622  }
1623  return -time_n;
1624  }
1625  return time_n;
1626 }
1627 
1628 
1629 
1630 //-----------------------------------------------------------
1631 // Constants
1632 //-----------------------------------------------------------
1633 
1634 const double DampFactor = 4;
1635 const int NDAMP = 10;
1636 
1637 
1638 //-----------------------------------------------------------
1639 // MultiNewton methods
1640 //-----------------------------------------------------------
1641 /**
1642  * L2 Norm of a delta in the solution
1643  *
1644  * The second argument has a default of false. However,
1645  * if true, then a table of the largest values is printed
1646  * out to standard output.
1647  */
1648 double BEulerInt::soln_error_norm(const double* const delta_y,
1649  bool printLargest)
1650 {
1651  int i;
1652  double sum_norm = 0.0, error;
1653  for (i = 0; i < m_neq; i++) {
1654  error = delta_y[i] / m_ewt[i];
1655  sum_norm += (error * error);
1656  }
1657  sum_norm = sqrt(sum_norm / m_neq);
1658  if (printLargest) {
1659  const int num_entries = 8;
1660  double dmax1, normContrib;
1661  int j;
1662  int* imax = mdp_alloc_int_1(num_entries, -1);
1663  printf("\t\tPrintout of Largest Contributors to norm "
1664  "of value (%g)\n", sum_norm);
1665  printf("\t\t I ysoln deltaY weightY "
1666  "Error_Norm**2\n");
1667  printf("\t\t ");
1668  print_line("-", 80);
1669  for (int jnum = 0; jnum < num_entries; jnum++) {
1670  dmax1 = -1.0;
1671  for (i = 0; i < m_neq; i++) {
1672  bool used = false;
1673  for (j = 0; j < jnum; j++) {
1674  if (imax[j] == i) {
1675  used = true;
1676  }
1677  }
1678  if (!used) {
1679  error = delta_y[i] / m_ewt[i];
1680  normContrib = sqrt(error * error);
1681  if (normContrib > dmax1) {
1682  imax[jnum] = i;
1683  dmax1 = normContrib;
1684  }
1685  }
1686  }
1687  i = imax[jnum];
1688  if (i >= 0) {
1689  printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e\n",
1690  i, m_y_n[i], delta_y[i], m_ewt[i], dmax1);
1691  }
1692  }
1693  printf("\t\t ");
1694  print_line("-", 80);
1695  mdp_safe_free((void**) &imax);
1696  }
1697  return sum_norm;
1698 }
1699 #ifdef DEBUG_HKM_JAC
1700 SquareMatrix jacBack();
1701 #endif
1702 /**************************************************************************
1703  *
1704  * doNewtonSolve():
1705  *
1706  * Compute the undamped Newton step. The residual function is
1707  * evaluated at the current time, t_n, at the current values of the
1708  * solution vector, m_y_n, and the solution time derivative, m_ydot_n,
1709  * but the Jacobian is not recomputed.
1710  */
1711 void BEulerInt::doNewtonSolve(double time_curr, double* y_curr,
1712  double* ydot_curr, double* delta_y,
1713  GeneralMatrix& jac, int loglevel)
1714 {
1715  int irow, jcol;
1716 
1717  m_func->evalResidNJ(time_curr, delta_t_n, y_curr,
1718  ydot_curr, delta_y, Base_ResidEval);
1719  m_nfe++;
1720  int sz = m_func->nEquations();
1721  for (int n = 0; n < sz; n++) {
1722  delta_y[n] = -delta_y[n];
1723  }
1724 
1725 
1726  /*
1727  * Column scaling -> We scale the columns of the Jacobian
1728  * by the nominal important change in the solution vector
1729  */
1730  if (m_colScaling) {
1731  if (!jac.factored()) {
1732  /*
1733  * Go get new scales
1734  */
1735  setColumnScales();
1736 
1737  /*
1738  * Scale the new Jacobian
1739  */
1740  double* jptr = &(*(jac.begin()));
1741  for (jcol = 0; jcol < m_neq; jcol++) {
1742  for (irow = 0; irow < m_neq; irow++) {
1743  *jptr *= m_colScales[jcol];
1744  jptr++;
1745  }
1746  }
1747  }
1748  }
1749 
1750  if (m_matrixConditioning) {
1751  if (jac.factored()) {
1752  m_func->matrixConditioning(0, sz, delta_y);
1753  } else {
1754  double* jptr = &(*(jac.begin()));
1755  m_func->matrixConditioning(jptr, sz, delta_y);
1756  }
1757  }
1758 
1759  /*
1760  * row sum scaling -> Note, this is an unequivocal success
1761  * at keeping the small numbers well balanced and
1762  * nonnegative.
1763  */
1764  if (m_rowScaling) {
1765  if (! jac.factored()) {
1766  /*
1767  * Ok, this is ugly. jac.begin() returns an vector<double> iterator
1768  * to the first data location.
1769  * Then &(*()) reverts it to a double *.
1770  */
1771  double* jptr = &(*(jac.begin()));
1772  for (irow = 0; irow < m_neq; irow++) {
1773  m_rowScales[irow] = 0.0;
1774  }
1775  for (jcol = 0; jcol < m_neq; jcol++) {
1776  for (irow = 0; irow < m_neq; irow++) {
1777  m_rowScales[irow] += fabs(*jptr);
1778  jptr++;
1779  }
1780  }
1781 
1782  jptr = &(*(jac.begin()));
1783  for (jcol = 0; jcol < m_neq; jcol++) {
1784  for (irow = 0; irow < m_neq; irow++) {
1785  *jptr /= m_rowScales[irow];
1786  jptr++;
1787  }
1788  }
1789  }
1790  for (irow = 0; irow < m_neq; irow++) {
1791  delta_y[irow] /= m_rowScales[irow];
1792  }
1793  }
1794 
1795 #ifdef DEBUG_HKM_JAC
1796  bool printJacContributions = false;
1797  if (m_time_step_num > 304) {
1798  printJacContributions = false;
1799  }
1800  int focusRow = 10;
1801  int numRows = 2;
1802  double RRow[2];
1803  bool freshJac = true;
1804  RRow[0] = delta_y[focusRow];
1805  RRow[1] = delta_y[focusRow+1];
1806  double Pcutoff = 1.0E-70;
1807  if (!jac.m_factored) {
1808  jacBack = jac;
1809  } else {
1810  freshJac = false;
1811  }
1812 #endif
1813  /*
1814  * Solve the system -> This also involves inverting the
1815  * matrix
1816  */
1817  (void) jac.solve(delta_y);
1818 
1819 
1820  /*
1821  * reverse the column scaling if there was any.
1822  */
1823  if (m_colScaling) {
1824  for (irow = 0; irow < m_neq; irow++) {
1825  delta_y[irow] *= m_colScales[irow];
1826  }
1827  }
1828 
1829 #ifdef DEBUG_HKM_JAC
1830  if (printJacContributions) {
1831  for (int iNum = 0; iNum < numRows; iNum++) {
1832  if (iNum > 0) {
1833  focusRow++;
1834  }
1835  double dsum = 0.0;
1836  vector_fp& Jdata = jacBack.data();
1837  double dRow = Jdata[m_neq * focusRow + focusRow];
1838  printf("\n Details on delta_Y for row %d \n", focusRow);
1839  printf(" Value before = %15.5e, delta = %15.5e,"
1840  "value after = %15.5e\n", y_curr[focusRow],
1841  delta_y[focusRow],
1842  y_curr[focusRow] + delta_y[focusRow]);
1843  if (!freshJac) {
1844  printf(" Old Jacobian\n");
1845  }
1846  printf(" col delta_y aij "
1847  "contrib \n");
1848  printf("--------------------------------------------------"
1849  "---------------------------------------------\n");
1850  printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
1851  focusRow, delta_y[focusRow],
1852  dRow, RRow[iNum] / dRow, RRow[iNum]);
1853  dsum += RRow[iNum] / dRow;
1854  for (int ii = 0; ii < m_neq; ii++) {
1855  if (ii != focusRow) {
1856  double aij = Jdata[m_neq * ii + focusRow];
1857  double contrib = aij * delta_y[ii] * (-1.0) / dRow;
1858  dsum += contrib;
1859  if (fabs(contrib) > Pcutoff) {
1860  printf("%6d %15.5e %15.5e %15.5e\n", ii,
1861  delta_y[ii] , aij, contrib);
1862  }
1863  }
1864  }
1865  printf("--------------------------------------------------"
1866  "---------------------------------------------\n");
1867  printf(" %15.5e %15.5e\n",
1868  delta_y[focusRow], dsum);
1869  }
1870  }
1871 
1872 #endif
1873 
1875 }
1876 
1877 //================================================================================================
1878 // Bound the Newton step while relaxing the solution
1879 /*
1880  * Return the factor by which the undamped Newton step 'step0'
1881  * must be multiplied in order to keep all solution components in
1882  * all domains between their specified lower and upper bounds.
1883  * Other bounds may be applied here as well.
1884  *
1885  * Currently the bounds are hard coded into this routine:
1886  *
1887  * Minimum value for all variables: - 0.01 * m_ewt[i]
1888  * Maximum value = none.
1889  *
1890  * Thus, this means that all solution components are expected
1891  * to be numerical greater than zero in the limit of time step
1892  * truncation errors going to zero.
1893  *
1894  * Delta bounds: The idea behind these is that the Jacobian
1895  * couldn't possibly be representative if the
1896  * variable is changed by a lot. (true for
1897  * nonlinear systems, false for linear systems)
1898  * Maximum increase in variable in any one newton iteration:
1899  * factor of 2
1900  * Maximum decrease in variable in any one newton iteration:
1901  * factor of 5
1902  *
1903  * @param y Current value of the solution
1904  * @param step0 Current raw step change in y[]
1905  * @param loglevel Log level. This routine produces output if loglevel
1906  * is greater than one
1907  *
1908  * @return Returns the damping coefficient
1909  */
1910 double BEulerInt::boundStep(const double* const y,
1911  const double* const step0, int loglevel)
1912 {
1913  int i, i_lower = -1, ifbd = 0, i_fbd = 0;
1914  double fbound = 1.0, f_lowbounds = 1.0, f_delta_bounds = 1.0;
1915  double ff, y_new, ff_alt;
1916  for (i = 0; i < m_neq; i++) {
1917  y_new = y[i] + step0[i];
1918  if ((y_new < (-0.01 * m_ewt[i])) && y[i] >= 0.0) {
1919  ff = 0.9 * (y[i] / (y[i] - y_new));
1920  if (ff < f_lowbounds) {
1921  f_lowbounds = ff;
1922  i_lower = i;
1923  }
1924  }
1925  /*
1926  * Now do a delta bounds
1927  * Increase variables by a factor of 2 only
1928  * decrease variables by a factor of 5 only
1929  */
1930  ff = 1.0;
1931  if ((fabs(y_new) > 2.0 * fabs(y[i])) &&
1932  (fabs(y_new-y[i]) > m_ewt[i])) {
1933  ff = fabs(y[i]/(y_new - y[i]));
1934  ff_alt = fabs(m_ewt[i] / (y_new - y[i]));
1935  ff = std::max(ff, ff_alt);
1936  ifbd = 1;
1937  }
1938  if ((fabs(5.0 * y_new) < fabs(y[i])) &&
1939  (fabs(y_new - y[i]) > m_ewt[i])) {
1940  ff = y[i]/(y_new-y[i]) * (1.0 - 5.0)/5.0;
1941  ff_alt = fabs(m_ewt[i] / (y_new - y[i]));
1942  ff = std::max(ff, ff_alt);
1943  ifbd = 0;
1944  }
1945  if (ff < f_delta_bounds) {
1946  f_delta_bounds = ff;
1947  i_fbd = ifbd;
1948  }
1949  f_delta_bounds = std::min(f_delta_bounds, ff);
1950  }
1951  fbound = std::min(f_lowbounds, f_delta_bounds);
1952  /*
1953  * Report on any corrections
1954  */
1955  if (loglevel > 1) {
1956  if (fbound != 1.0) {
1957  if (f_lowbounds < f_delta_bounds) {
1958  printf("\t\tboundStep: Variable %d causing lower bounds "
1959  "damping of %g\n",
1960  i_lower, f_lowbounds);
1961  } else {
1962  if (ifbd) {
1963  printf("\t\tboundStep: Decrease of Variable %d causing "
1964  "delta damping of %g\n",
1965  i_fbd, f_delta_bounds);
1966  } else {
1967  printf("\t\tboundStep: Increase of variable %d causing"
1968  "delta damping of %g\n",
1969  i_fbd, f_delta_bounds);
1970  }
1971  }
1972  }
1973  }
1974  return fbound;
1975 }
1976 //================================================================================================
1977 /**************************************************************************
1978  *
1979  * dampStep():
1980  *
1981  * On entry, step0 must contain an undamped Newton step for the
1982  * solution x0. This method attempts to find a damping coefficient
1983  * such that the next undamped step would have a norm smaller than
1984  * that of step0. If successful, the new solution after taking the
1985  * damped step is returned in y1, and the undamped step at y1 is
1986  * returned in step1.
1987  */
1988 int BEulerInt::dampStep(double time_curr, const double* y0,
1989  const double* ydot0, const double* step0,
1990  double* y1, double* ydot1, double* step1,
1991  double& s1, GeneralMatrix& jac,
1992  int& loglevel, bool writetitle,
1993  int& num_backtracks)
1994 {
1995 
1996 
1997  // Compute the weighted norm of the undamped step size step0
1998  double s0 = soln_error_norm(step0);
1999 
2000  // Compute the multiplier to keep all components in bounds
2001  // A value of one indicates that there is no limitation
2002  // on the current step size in the nonlinear method due to
2003  // bounds constraints (either negative values of delta
2004  // bounds constraints.
2005  double fbound = boundStep(y0, step0, loglevel);
2006 
2007  // if fbound is very small, then y0 is already close to the
2008  // boundary and step0 points out of the allowed domain. In
2009  // this case, the Newton algorithm fails, so return an error
2010  // condition.
2011  if (fbound < 1.e-10) {
2012  if (loglevel > 1) {
2013  printf("\t\t\tdampStep: At limits.\n");
2014  }
2015  return -3;
2016  }
2017 
2018  //--------------------------------------------
2019  // Attempt damped step
2020  //--------------------------------------------
2021 
2022  // damping coefficient starts at 1.0
2023  double damp = 1.0;
2024  int j, m;
2025  double ff;
2026  num_backtracks = 0;
2027  for (m = 0; m < NDAMP; m++) {
2028 
2029  ff = fbound*damp;
2030 
2031  // step the solution by the damped step size
2032  /*
2033  * Whenever we update the solution, we must also always
2034  * update the time derivative.
2035  */
2036  for (j = 0; j < m_neq; j++) {
2037  y1[j] = y0[j] + ff*step0[j];
2038  // HKM setting intermediate y's to zero was a tossup.
2039  // slightly different, equivalent results
2040  //#ifdef DEBUG_HKM
2041  // y1[j] = MAX(0.0, y1[j]);
2042  //#endif
2043  }
2044  calc_ydot(m_order, y1, ydot1);
2045 
2046  // compute the next undamped step, step1[], that would result
2047  // if y1[] were accepted.
2048 
2049  doNewtonSolve(time_curr, y1, ydot1, step1, jac, loglevel);
2050 
2051 #ifdef DEBUG_HKM
2052  for (j = 0; j < m_neq; j++) {
2053  checkFinite(step1[j]);
2054  checkFinite(y1[j]);
2055  }
2056 #endif
2057  // compute the weighted norm of step1
2058  s1 = soln_error_norm(step1);
2059 
2060  // write log information
2061  if (loglevel > 3) {
2062  print_solnDelta_norm_contrib((const double*) step0,
2063  "DeltaSolnTrial",
2064  (const double*) step1,
2065  "DeltaSolnTrialTest",
2066  "dampNewt: Important Entries for "
2067  "Weighted Soln Updates:",
2068  y0, y1, ff, 5);
2069  }
2070  if (loglevel > 1) {
2071  printf("\t\t\tdampNewt: s0 = %g, s1 = %g, fbound = %g,"
2072  "damp = %g\n", s0, s1, fbound, damp);
2073  }
2074 #ifdef DEBUG_HKM
2075  if (loglevel > 2) {
2076  if (s1 > 1.00000001 * s0 && s1 > 1.0E-5) {
2077  printf("WARNING: Possible Jacobian Problem "
2078  "-> turning on more debugging for this step!!!\n");
2079  print_solnDelta_norm_contrib((const double*) step0,
2080  "DeltaSolnTrial",
2081  (const double*) step1,
2082  "DeltaSolnTrialTest",
2083  "dampNewt: Important Entries for "
2084  "Weighted Soln Updates:",
2085  y0, y1, ff, 5);
2086  loglevel = 4;
2087  }
2088  }
2089 #endif
2090 
2091  // if the norm of s1 is less than the norm of s0, then
2092  // accept this damping coefficient. Also accept it if this
2093  // step would result in a converged solution. Otherwise,
2094  // decrease the damping coefficient and try again.
2095 
2096  if (s1 < 1.0E-5 || s1 < s0) {
2097  if (loglevel > 2) {
2098  if (s1 > s0) {
2099  if (s1 > 1.0) {
2100  printf("\t\t\tdampStep: current trial step and damping"
2101  " coefficient accepted because test step < 1\n");
2102  printf("\t\t\t s1 = %g, s0 = %g\n", s1, s0);
2103  }
2104  }
2105  }
2106  break;
2107  } else {
2108  if (loglevel > 1) {
2109  printf("\t\t\tdampStep: current step rejected: (s1 = %g > "
2110  "s0 = %g)", s1, s0);
2111  if (m < (NDAMP-1)) {
2112  printf(" Decreasing damping factor and retrying");
2113  } else {
2114  printf(" Giving up!!!");
2115  }
2116  printf("\n");
2117  }
2118  }
2119  num_backtracks++;
2120  damp /= DampFactor;
2121  }
2122 
2123  // If a damping coefficient was found, return 1 if the
2124  // solution after stepping by the damped step would represent
2125  // a converged solution, and return 0 otherwise. If no damping
2126  // coefficient could be found, return -2.
2127  if (m < NDAMP) {
2128  if (s1 > 1.0) {
2129  return 0;
2130  } else {
2131  return 1;
2132  }
2133  } else {
2134  if (s1 < 0.5 && (s0 < 0.5)) {
2135  return 1;
2136  }
2137  if (s1 < 1.0) {
2138  return 0;
2139  }
2140  return -2;
2141  }
2142 }
2143 //================================================================================================
2144 // Solve a nonlinear system
2145 /*
2146  * Find the solution to F(X, xprime) = 0 by damped Newton iteration. On
2147  * entry, y_comm[] contains an initial estimate of the solution and
2148  * ydot_comm[] contains an estimate of the derivative.
2149  * On successful return, y_comm[] contains the converged solution
2150  * and ydot_comm[] contains the derivative
2151  *
2152  *
2153  * @param y_comm[] Contains the input solution. On output y_comm[] contains
2154  * the converged solution
2155  * @param ydot_comm Contains the input derivative solution. On output y_comm[] contains
2156  * the converged derivative solution
2157  * @param CJ Inverse of the time step
2158  * @param time_curr Current value of the time
2159  * @param jac Jacobian
2160  * @param num_newt_its number of newton iterations
2161  * @param num_linear_solves number of linear solves
2162  * @param num_backtracks number of backtracs
2163  * @param loglevel Log level
2164  */
2165 int BEulerInt::solve_nonlinear_problem(double* const y_comm,
2166  double* const ydot_comm, double CJ,
2167  double time_curr,
2168  GeneralMatrix& jac,
2169  int& num_newt_its,
2170  int& num_linear_solves,
2171  int& num_backtracks,
2172  int loglevel)
2173 {
2174  int m = 0;
2175  bool forceNewJac = false;
2176  double s1=1.e30;
2177 
2178  double* y_curr = mdp_alloc_dbl_1(m_neq, 0.0);
2179  double* ydot_curr = mdp_alloc_dbl_1(m_neq, 0.0);
2180  double* stp = mdp_alloc_dbl_1(m_neq, 0.0);
2181  double* stp1 = mdp_alloc_dbl_1(m_neq, 0.0);
2182  double* y_new = mdp_alloc_dbl_1(m_neq, 0.0);
2183  double* ydot_new = mdp_alloc_dbl_1(m_neq, 0.0);
2184 
2185  mdp_copy_dbl_1(y_curr, y_comm, m_neq);
2186  mdp_copy_dbl_1(ydot_curr, ydot_comm, m_neq);
2187 
2188  bool frst = true;
2189  num_newt_its = 0;
2190  num_linear_solves = - m_numTotalLinearSolves;
2191  num_backtracks = 0;
2192  int i_backtracks;
2193 
2194  while (1 > 0) {
2195 
2196  /*
2197  * Increment Newton Solve counter
2198  */
2200  num_newt_its++;
2201 
2202 
2203  if (loglevel > 1) {
2204  printf("\t\tSolve_Nonlinear_Problem: iteration %d:\n",
2205  num_newt_its);
2206  }
2207 
2208  // Check whether the Jacobian should be re-evaluated.
2209 
2210  forceNewJac = true;
2211 
2212  if (forceNewJac) {
2213  if (loglevel > 1) {
2214  printf("\t\t\tGetting a new Jacobian and solving system\n");
2215  }
2216  beuler_jac(jac, m_resid, time_curr, CJ, y_curr, ydot_curr,
2217  num_newt_its);
2218  } else {
2219  if (loglevel > 1) {
2220  printf("\t\t\tSolving system with old jacobian\n");
2221  }
2222  }
2223 
2224  // compute the undamped Newton step
2225  doNewtonSolve(time_curr, y_curr, ydot_curr, stp, jac, loglevel);
2226 
2227  // damp the Newton step
2228  m = dampStep(time_curr, y_curr, ydot_curr, stp, y_new, ydot_new,
2229  stp1, s1, jac, loglevel, frst, i_backtracks);
2230  frst = false;
2231  num_backtracks += i_backtracks;
2232 
2233  /*
2234  * Impose the minimum number of newton iterations critera
2235  */
2236  if (num_newt_its < m_min_newt_its) {
2237  if (m == 1) {
2238  m = 0;
2239  }
2240  }
2241  /*
2242  * Impose max newton iteration
2243  */
2244  if (num_newt_its > 20) {
2245  m = -1;
2246  if (loglevel > 1) {
2247  printf("\t\t\tDampnewton unsuccessful (max newts exceeded) sfinal = %g\n", s1);
2248  }
2249  }
2250 
2251  if (loglevel > 1) {
2252  if (m == 1) {
2253  printf("\t\t\tDampNewton iteration successful, nonlin "
2254  "converged sfinal = %g\n", s1);
2255  } else if (m == 0) {
2256  printf("\t\t\tDampNewton iteration successful, get new"
2257  "direction, sfinal = %g\n", s1);
2258  } else {
2259  printf("\t\t\tDampnewton unsuccessful sfinal = %g\n", s1);
2260  }
2261  }
2262 
2263  // If we are converged, then let's use the best solution possible
2264  // for an end result. We did a resolve in dampStep(). Let's update
2265  // the solution to reflect that.
2266  // HKM 5/16 -> Took this out, since if the last step was a
2267  // damped step, then adding stp1[j] is undamped, and
2268  // may lead to oscillations. It kind of defeats the
2269  // purpose of dampStep() anyway.
2270  // if (m == 1) {
2271  // for (int j = 0; j < m_neq; j++) {
2272  // y_new[j] += stp1[j];
2273  // HKM setting intermediate y's to zero was a tossup.
2274  // slightly different, equivalent results
2275  // #ifdef DEBUG_HKM
2276  // y_new[j] = MAX(0.0, y_new[j]);
2277  // #endif
2278  // }
2279  // }
2280 
2281  bool m_filterIntermediate = false;
2282  if (m_filterIntermediate) {
2283  if (m == 0) {
2284  (void) filterNewStep(time_n, y_new, ydot_new);
2285  }
2286  }
2287  // Exchange new for curr solutions
2288  if (m == 0 || m == 1) {
2289  mdp_copy_dbl_1(y_curr, y_new, m_neq);
2290  calc_ydot(m_order, y_curr, ydot_curr);
2291  }
2292 
2293  // convergence
2294  if (m == 1) {
2295  goto done;
2296  }
2297 
2298  // If dampStep fails, first try a new Jacobian if an old
2299  // one was being used. If it was a new Jacobian, then
2300  // return -1 to signify failure.
2301  else if (m < 0) {
2302  goto done;
2303  }
2304  }
2305 
2306 done:
2307  // Copy into the return vectors
2308  mdp_copy_dbl_1(y_comm, y_curr, m_neq);
2309  mdp_copy_dbl_1(ydot_comm, ydot_curr, m_neq);
2310  // Increment counters
2311  num_linear_solves += m_numTotalLinearSolves;
2312  // Free memory
2313  mdp_safe_free((void**) &y_curr);
2314  mdp_safe_free((void**) &ydot_curr);
2315  mdp_safe_free((void**) &stp);
2316  mdp_safe_free((void**) &stp1);
2317  mdp_safe_free((void**) &y_new);
2318  mdp_safe_free((void**) &ydot_new);
2319 
2320  double time_elapsed = 0.0;
2321  if (loglevel > 1) {
2322  if (m == 1) {
2323  printf("\t\tNonlinear problem solved successfully in "
2324  "%d its, time elapsed = %g sec\n",
2325  num_newt_its, time_elapsed);
2326  }
2327  }
2328  return m;
2329 }
2330 //================================================================================================
2331 /*
2332  *
2333  *
2334  */
2335 void BEulerInt::
2336 print_solnDelta_norm_contrib(const double* const solnDelta0,
2337  const char* const s0,
2338  const double* const solnDelta1,
2339  const char* const s1,
2340  const char* const title,
2341  const double* const y0,
2342  const double* const y1,
2343  double damp,
2344  int num_entries)
2345 {
2346  int i, j, jnum;
2347  bool used;
2348  double dmax0, dmax1, error, rel_norm;
2349  printf("\t\t%s currentDamp = %g\n", title, damp);
2350  printf("\t\t I ysoln %10s ysolnTrial "
2351  "%10s weight relSoln0 relSoln1\n", s0, s1);
2352  int* imax = mdp_alloc_int_1(num_entries, -1);
2353  printf("\t\t ");
2354  print_line("-", 90);
2355  for (jnum = 0; jnum < num_entries; jnum++) {
2356  dmax1 = -1.0;
2357  for (i = 0; i < m_neq; i++) {
2358  used = false;
2359  for (j = 0; j < jnum; j++) {
2360  if (imax[j] == i) {
2361  used = true;
2362  }
2363  }
2364  if (!used) {
2365  error = solnDelta0[i] / m_ewt[i];
2366  rel_norm = sqrt(error * error);
2367  error = solnDelta1[i] / m_ewt[i];
2368  rel_norm += sqrt(error * error);
2369  if (rel_norm > dmax1) {
2370  imax[jnum] = i;
2371  dmax1 = rel_norm;
2372  }
2373  }
2374  }
2375  if (imax[jnum] >= 0) {
2376  i = imax[jnum];
2377  error = solnDelta0[i] / m_ewt[i];
2378  dmax0 = sqrt(error * error);
2379  error = solnDelta1[i] / m_ewt[i];
2380  dmax1 = sqrt(error * error);
2381  printf("\t\t %4d %12.4e %12.4e %12.4e %12.4e "
2382  "%12.4e %12.4e %12.4e\n",
2383  i, y0[i], solnDelta0[i], y1[i],
2384  solnDelta1[i], m_ewt[i], dmax0, dmax1);
2385  }
2386  }
2387  printf("\t\t ");
2388  print_line("-", 90);
2389  mdp_safe_free((void**) &imax);
2390 }
2391 //===============================================================================================
2392 
2393 } // End of namespace Cantera
2394