Cantera  2.0
BEulerInt.h
Go to the documentation of this file.
1 /**
2  * @file BEulerInt.h
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 #ifndef CT_BEULERINT_H
12 #define CT_BEULERINT_H
13 
14 #include "cantera/base/ct_defs.h"
15 #include "cantera/numerics/ctlapack.h"
16 #include "cantera/base/utilities.h"
17 #include "cantera/base/ctexceptions.h"
18 
19 
20 #include "cantera/numerics/Integrator.h"
21 #include "cantera/numerics/ResidJacEval.h"
22 
23 #include "cantera/numerics/GeneralMatrix.h"
24 #include "cantera/numerics/NonlinearSolver.h"
25 
26 #include "cantera/base/mdp_allo.h"
27 
28 #define OPT_SIZE 10
29 
30 #define SUCCESS 0
31 #define FAILURE 1
32 
33 #define STEADY 0
34 #define TRANSIENT 1
35 
36 namespace Cantera
37 {
38 
39 enum BEulerMethodType {
40  BEulerFixedStep,
41  BEulerVarStep
42 };
43 
44 /**
45  * Exception class thrown when a BEuler error is encountered.
46  */
47 class BEulerErr : public CanteraError
48 {
49 public:
50  BEulerErr(std::string msg);
51 };
52 
53 
54 #define BEULER_JAC_ANAL 2
55 #define BEULER_JAC_NUM 1
56 
57 /**
58  * Wrapper class for 'beuler' integrator
59  * We derive the class from the class Integrator
60  */
61 class BEulerInt : public Integrator
62 {
63 
64 public:
65 
66  //! The default constructor doesn't take an argument.
67  BEulerInt();
68  //! Destructor
69  virtual ~BEulerInt();
70  virtual void setTolerances(double reltol, size_t n, double* abstol);
71  virtual void setTolerances(double reltol, double abstol);
72  virtual void setProblemType(int probtype);
73  virtual void initializeRJE(double t0, ResidJacEval& func);
74  virtual void reinitializeRJE(double t0, ResidJacEval& func);
75  virtual double integrateRJE(double tout, double tinit = 0.0);
76  virtual doublereal step(double tout);
77  virtual void setSolnWeights();
78  virtual double& solution(size_t k) {
79  return m_y_n[k];
80  }
81  double* solution() {
82  return m_y_n;
83  }
84  int nEquations() const {
85  return m_neq;
86  }
87  virtual int nEvals() const;
88  virtual void setMethodBEMT(BEulerMethodType t);
89  virtual void setIterator(IterType t);
90  virtual void setMaxStep(double hmax);
91  virtual void setMaxNumTimeSteps(int);
92  virtual void setNumInitialConstantDeltaTSteps(int);
93 
94  void print_solnDelta_norm_contrib(const double* const soln0,
95  const char* const s0,
96  const double* const soln1,
97  const char* const s1,
98  const char* const title,
99  const double* const y0,
100  const double* const y1,
101  double damp,
102  int num_entries);
103 
104  virtual void setPrintSolnOptions(int printSolnStepInterval,
105  int printSolnNumberToTout,
106  int printSolnFirstSteps = 0,
107  bool dumpJacobians = false);
108  void setNonLinOptions(int min_newt_its = 0,
109  bool matrixConditioning = false,
110  bool colScaling = false,
111  bool rowScaling = true);
112  virtual void setPrintFlag(int print_flag);
113  virtual void setColumnScales();
114  /**
115  * calculate the solution error norm
116  */
117  virtual double soln_error_norm(const double* const,
118  bool printLargest = false);
119  virtual void setInitialTimeStep(double delta_t);
120 
121  void beuler_jac(GeneralMatrix&, double* const,
122  double, double, double* const, double* const, int);
123 
124 
125 protected:
126 
127  //! Internal routine that sets up the fixed length storage based on
128  //! the size of the problem to solve.
129  void internalMalloc();
130  /**
131  * Internal function to calculate the predicted solution
132  * at a time step.
133  */
134  void calc_y_pred(int);
135  /**
136  * Internal function to calculate the time derivative at the
137  * new step
138  */
139  void calc_ydot(int, double*, double*);
140  /**
141  * Internal function to calculate the time step truncation
142  * error for a predictor corrector time step
143  */
144  double time_error_norm();
145  /**
146  * Internal function to calculate the time step for the
147  * next step based on the time-truncation error on the
148  * current time step
149  */
150  double time_step_control(int m_order, double time_error_factor);
151 
152  //! Solve a nonlinear system
153  /*!
154  *
155  * Find the solution to F(X, xprime) = 0 by damped Newton iteration. On
156  * entry, y_comm[] contains an initial estimate of the solution and
157  * ydot_comm[] contains an estimate of the derivative.
158  * On successful return, y_comm[] contains the converged solution
159  * and ydot_comm[] contains the derivative
160  *
161  *
162  * @param y_comm[] Contains the input solution. On output y_comm[] contains
163  * the converged solution
164  * @param ydot_comm Contains the input derivative solution. On output y_comm[] contains
165  * the converged derivative solution
166  * @param CJ Inverse of the time step
167  * @param time_curr Current value of the time
168  * @param jac Jacobian
169  * @param num_newt_its number of newton iterations
170  * @param num_linear_solves number of linear solves
171  * @param num_backtracks number of backtracs
172  * @param loglevel Log level
173  */
174  int solve_nonlinear_problem(double* const y_comm,
175  double* const ydot_comm, double CJ,
176  double time_curr,
177  GeneralMatrix& jac,
178  int& num_newt_its,
179  int& num_linear_solves,
180  int& num_backtracks,
181  int loglevel);
182 
183  /**
184  * Compute the undamped Newton step. The residual function is
185  * evaluated at x, but the Jacobian is not recomputed.
186  */
187  void doNewtonSolve(double, double*, double*, double*,
188  GeneralMatrix&, int);
189 
190 
191  //! Bound the Newton step while relaxing the solution
192  /*!
193  * Return the factor by which the undamped Newton step 'step0'
194  * must be multiplied in order to keep all solution components in
195  * all domains between their specified lower and upper bounds.
196  * Other bounds may be applied here as well.
197  *
198  * Currently the bounds are hard coded into this routine:
199  *
200  * Minimum value for all variables: - 0.01 * m_ewt[i]
201  * Maximum value = none.
202  *
203  * Thus, this means that all solution components are expected
204  * to be numerical greater than zero in the limit of time step
205  * truncation errors going to zero.
206  *
207  * Delta bounds: The idea behind these is that the Jacobian
208  * couldn't possibly be representative if the
209  * variable is changed by a lot. (true for
210  * nonlinear systems, false for linear systems)
211  * Maximum increase in variable in any one newton iteration:
212  * factor of 2
213  * Maximum decrease in variable in any one newton iteration:
214  * factor of 5
215  *
216  * @param y Current value of the solution
217  * @param step0 Current raw step change in y[]
218  * @param loglevel Log level. This routine produces output if loglevel
219  * is greater than one
220  *
221  * @return Returns the damping coefficient
222  */
223  double boundStep(const double* const y, const double* const step0, int loglevel);
224 
225  /*
226  * Damp step
227  */
228  int dampStep(double, const double*, const double*,
229  const double*, double*, double*,
230  double*, double&, GeneralMatrix&, int&, bool, int&);
231 
232  /*
233  * Compute Residual Weights
234  */
235  void computeResidWts(GeneralMatrix& jac);
236 
237  /*
238  * Filter a new step
239  */
240  double filterNewStep(double, double*, double*);
241 
242  /*
243  * get the next time to print out
244  */
245  double getPrintTime(double time_current);
246 
247  /********************** Member data ***************************/
248  /*********************
249  * METHOD FLAGS
250  *********************/
251 
252  //! IterType is used to specify how the nonlinear equations are
253  //! to be relaxed at each time step.
254  IterType m_iter;
255  /**
256  * MethodType is used to specify how the time step is to be
257  * chosen. Currently, there are two choices, one is a fixed
258  * step method while the other is based on a predictor-corrector
259  * algorithm and a time-step truncation error tolerance.
260  */
261  BEulerMethodType m_method;
262  /**
263  * m_jacFormMethod determines how a matrix is formed.
264  */
265  int m_jacFormMethod;
266  /**
267  * m_rowScaling is a boolean. If true then row sum scaling
268  * of the Jacobian matrix is carried out when solving the
269  * linear systems.
270  */
271  bool m_rowScaling;
272  /**
273  * m_colScaling is a boolean. If true, then column scaling
274  * is performed on each solution of the linear system.
275  */
276  bool m_colScaling;
277  /**
278  * m_matrixConditioning is a boolean. If true, then the
279  * Jacobian and every rhs is multiplied by the inverse
280  * of a matrix that is suppose to reduce the condition
281  * number of the matrix. This is done before row scaling.
282  */
283  bool m_matrixConditioning;
284  /**
285  * If m_itol =1 then each component has an individual
286  * value of atol. If m_itol = 0, the all atols are equal.
287  */
288  int m_itol;
289  /**
290  * Relative time truncation error tolerances
291  */
292  double m_reltol;
293  /**
294  * Absolute time truncation error tolerances, when uniform
295  * for all variables.
296  */
297  double m_abstols;
298  /**
299  * Vector of absolute time truncation error tolerance
300  * when not uniform for all variables.
301  */
302  double* m_abstol;
303  /**
304  * Error Weights. This is a surprisingly important quantity.
305  */
306  double* m_ewt;
307 
308  //! Maximum step size
309  double m_hmax;
310  /**
311  * Maximum integration order
312  */
313  int m_maxord;
314  /**
315  * Current integration order
316  */
317  int m_order;
318  /**
319  * Time step number
320  */
321  int m_time_step_num;
322  int m_time_step_attempts;
323  /**
324  * Max time steps allowed
325  */
326  int m_max_time_step_attempts;
327  /**
328  * Number of initial time steps to take where the
329  * time truncation error tolerances are not checked. Instead
330  * the delta T is uniform
331  */
332  int m_numInitialConstantDeltaTSteps;
333  /**
334  * Failure Counter -> keeps track of the number
335  * of consequetive failures
336  */
337  int m_failure_counter;
338  /**
339  * Minimum Number of Newton Iterations per nonlinear step
340  * default = 0
341  */
342  int m_min_newt_its;
343  /************************
344  * PRINTING OPTIONS
345  ************************/
346  /**
347  * Step Interval at which to print out the solution
348  * default = 1;
349  * If set to zero, there is no printout
350  */
351  int m_printSolnStepInterval;
352  /**
353  * Number of evenly spaced printouts of the solution
354  * If zero, there is no printout from this option
355  * default 1
356  * If set to zero there is no printout.
357  */
358  int m_printSolnNumberToTout;
359 
360  /**
361  * Number of initial steps that the solution is
362  * printed out.
363  * default = 0
364  */
365  int m_printSolnFirstSteps;
366 
367  /**
368  * Dump Jacobians to disk
369  * default false
370  */
371  bool m_dumpJacobians;
372 
373  /*********************
374  * INTERNAL SOLUTION VALUES
375  *********************/
376  /**
377  * Number of equations in the ode integrator
378  */
379  int m_neq;
380  double* m_y_n;
381  double* m_y_nm1;
382  double* m_y_pred_n;
383  double* m_ydot_n;
384  double* m_ydot_nm1;
385  /************************
386  * TIME VARIABLES
387  ************************/
388  /**
389  * Initial time at the start of the integration
390  */
391  double m_t0;
392  /**
393  * Final time
394  */
395  double m_time_final;
396  /**
397  *
398  */
399  double time_n;
400  double time_nm1;
401  double time_nm2;
402  double delta_t_n;
403  double delta_t_nm1;
404  double delta_t_nm2;
405  double delta_t_np1;
406  /**
407  * Maximum permissible time step
408  */
409  double delta_t_max;
410 
411 
412  double* m_resid;
413  double* m_residWts;
414  double* m_wksp;
415  ResidJacEval* m_func;
416  double* m_rowScales;
417  double* m_colScales;
418 
419  /**
420  * Pointer to the jacobian representing the
421  * time dependent problem
422  */
423  GeneralMatrix* tdjac_ptr;
424  /**
425  * Determines the level of printing for each time
426  * step.
427  * 0 -> absolutely nothing is printed for
428  * a single time step.
429  * 1 -> One line summary per time step
430  * 2 -> short description, points of interest
431  * 3 -> Lots printed per time step (default)
432  */
433  int m_print_flag;
434  /***************************************************************************
435  * COUNTERS OF VARIOUS KINDS
436  ***************************************************************************/
437  /**
438  * Number of function evaluations
439  */
440  int m_nfe;
441  /**
442  * Number of Jacobian Evaluations and
443  * factorization steps (they are the same)
444  */
445  int m_nJacEval;
446  /**
447  * Number of total newton iterations
448  */
449  int m_numTotalNewtIts;
450  /**
451  * Total number of linear iterations
452  */
453  int m_numTotalLinearSolves;
454  /**
455  * Total number of convergence failures.
456  */
457  int m_numTotalConvFails;
458  /**
459  * Total Number of time truncation error failures
460  */
461  int m_numTotalTruncFails;
462  /*
463  *
464  */
465  int num_failures;
466 };
467 
468 } // namespace
469 
470 #endif // CT_BEULER
Generated by   doxygen 1.8.2