Cantera  2.1.2
solveProb.h
Go to the documentation of this file.
1 /**
2  * @file solveProb.h Header file for implicit nonlinear solver with the option
3  * of a pseudotransient (see \ref numerics and class \link
5  */
6
7 /*
8  * Copyright 2004 Sandia Corporation. Under the terms of Contract
9  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
10  * retains certain rights in this software.
11  * See file License.txt for licensing information.
12  */
13
14 #ifndef SOLVEPROB_H
15 #define SOLVEPROB_H
16 /**
17  * @defgroup solverGroup Solvers for Equation Systems
18  */
19
20 #include "cantera/base/Array.h"
21 #include "ResidEval.h"
22
23 //! Solution Methods
24 /*!
25  * Flag to specify the solution method
26  *
27  * 1: SOLVEPROB_INITIALIZE = This assumes that the initial guess supplied to the
28  * routine is far from the correct one. Substantial
29  * work plus transient time-stepping is to be expected
30  * to find a solution.
31  * 2: SOLVEPROB_RESIDUAL = Need to solve the surface problem in order to
32  * calculate the surface fluxes of gas-phase species.
33  * (Can expect a moderate change in the solution
34  * vector -> try to solve the system by direct
35  * methods
36  * with no damping first -> then, try time-stepping
37  * if the first method fails)
38  * A "time_scale" supplied here is used in the
39  * algorithm to determine when to shut off
40  * time-stepping.
41  * 3: SOLVEPROB_JACOBIAN = Calculation of the surface problem is due to the
42  * need for a numerical jacobian for the gas-problem.
43  * The solution is expected to be very close to the
44  * initial guess, and accuracy is needed.
45  * 4: SOLVEPROB_TRANSIENT = The transient calculation is performed here for an
46  * amount of time specified by "time_scale". It is
47  * not guaranteed to be time-accurate - just stable
48  * and fairly fast. The solution after del_t time is
49  * returned, whether it's converged to a steady
50  * state or not.
51  */
52 const int SOLVEPROB_INITIALIZE = 1;
53 const int SOLVEPROB_RESIDUAL = 2;
54 const int SOLVEPROB_JACOBIAN = 3;
55 const int SOLVEPROB_TRANSIENT = 4;
56
57 namespace Cantera
58 {
59 //! Method to solve a pseudo steady state of a nonlinear problem
60 /*!
61  * The following class handles the solution of a nonlinear problem.
62  *
63  * Res_ss(C) = - Res(C) = 0
64  *
65  * Optionally a pseudo transient algorithm may be used to relax the residual if
66  * it is available.
67  *
68  * Res_td(C) = dC/dt - Res(C) = 0;
69  *
70  * Res_ss(C) is the steady state residual to be solved. Res_td(C) is the
71  * time dependent residual which leads to the steady state residual.
72  *
73  * Solution Method
74  *
75  * This routine is typically used within a residual calculation in a large code.
76  * It's typically invoked millions of times for large calculations, and it must
77  * work every time. Therefore, requirements demand that it be robust but also
78  * efficient.
79  *
80  * The solution methodology is largely determined by the <TT>ifunc</TT> parameter,
81  * that is input to the solution object. This parameter may have the following
82  * 4 values:
83  *
84  * 1: SOLVEPROB_INITIALIZE = This assumes that the initial guess supplied to the
85  * routine is far from the correct one. Substantial
86  * work plus transient time-stepping is to be expected
87  * to find a solution.
88  *
89  * 2: SOLVEPROB_RESIDUAL = Need to solve the nonlinear problem in order to
90  * calculate quantities for a residual calculation
91  * (Can expect a moderate change in the solution
92  * vector -> try to solve the system by direct methods
93  * with no damping first -> then, try time-stepping
94  * if the first method fails)
95  * A "time_scale" supplied here is used in the
96  * algorithm to determine when to shut off
97  * time-stepping.
98  *
99  * 3: SOLVEPROB_JACOBIAN = Calculation of the surface problem is due to the
100  * need for a numerical jacobian for the gas-problem.
101  * The solution is expected to be very close to the
102  * initial guess, and extra accuracy is needed because
103  * solution variables have been delta'd from
104  * nominal values to create jacobian entries.
105  *
106  * 4: SOLVEPROB_TRANSIENT = The transient calculation is performed here for an
107  * amount of time specified by "time_scale". It is
108  * not guaranteed to be time-accurate - just stable
109  * and fairly fast. The solution after del_t time is
110  * returned, whether it's converged to a steady
111  * state or not. This is a poor man's time stepping
112  * algorithm.
113  *
114  * Pseudo time stepping algorithm:
115  * The time step is determined from sdot[], so that the time step
116  * doesn't ever change the value of a variable by more than 100%.
117  *
118  * This algorithm does use a damped Newton's method to relax the equations.
119  * Damping is based on a "delta damping" technique. The solution unknowns
120  * are not allowed to vary too much between iterations.
121  *
122  * EXTRA_ACCURACY:A constant that is the ratio of the required update norm in
123  * this Newton iteration compared to that in the nonlinear solver.
124  * A value of 0.1 is used so surface species are safely overconverged.
125  *
126  * Functions called:
127  *----------------------------------------------------------------------------
128  *
129  * ct_dgetrf -- First half of LAPACK direct solve of a full Matrix
130  *
131  * ct_dgetrs -- Second half of LAPACK direct solve of a full matrix. Returns
132  * solution vector in the right-hand-side vector, resid.
133  *
134  *----------------------------------------------------------------------------
135  *
136  * @ingroup solverGroup
137  */
139 {
140 public:
141
142  //! Constructor for the object
143  solveProb(ResidEval* resid);
144
145  virtual ~solveProb();
146
147 private:
148
149  //! Unimplemented private copy constructor
150  solveProb(const solveProb& right);
151
152  //! Unimplemented private assignment operator
153  solveProb& operator=(const solveProb& right);
154
155 public:
156  //! Main routine that actually calculates the pseudo steady state
157  //! of the surface problem
158  /*!
159  * The actual converged solution is returned as part of the
160  * internal state of the InterfaceKinetics objects.
161  *
162  * @param ifunc Determines the type of solution algorithm to be
163  * used. Possible values are SOLVEPROB_INITIALIZE ,
164  * SOLVEPROB_RESIDUAL SOLVEPROB_JACOBIAN SOLVEPROB_TRANSIENT .
165  *
166  * @param time_scale Time over which to integrate the surface equations,
167  * where applicable
168  *
169  * @param reltol Relative tolerance to use
170  *
171  * @return Returns 1 if the surface problem is successfully solved.
172  * Returns -1 if the surface problem wasn't solved successfully.
173  * Note the actual converged solution is returned as part of the
174  * internal state of the InterfaceKinetics objects.
175  */
176  int solve(int ifunc, doublereal time_scale, doublereal reltol);
177
178  //! Report the current state of the solution
179  /*!
180  * @param[out] CSoln solution vector for the nonlinear problem
181  */
182  virtual void reportState(doublereal* const CSoln) const;
183
184  //! Set the bottom and top bounds on the solution vector
185  /*!
186  * The default is for the bottom is 0.0, while the default for the top is 1.0
187  *
188  * @param botBounds Vector of bottom bounds
189  * @param topBounds vector of top bounds
190  */
191  virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]);
192
193  void setAtol(const doublereal atol[]);
194  void setAtolConst(const doublereal atolconst);
195
196 private:
197  //! Printing routine that gets called at the start of every invocation
198  virtual void print_header(int ioflag, int ifunc, doublereal time_scale,
199  doublereal reltol,
200  doublereal netProdRate[]);
201
202 #ifdef DEBUG_SOLVEPROB
203  //! Prints out the residual and Jacobian
204  virtual void printResJac(int ioflag, int neq, const Array2D& Jac,
205  doublereal resid[], doublereal wtResid[], doublereal norm);
206 #endif
207
208  //! Printing routine that gets called after every iteration
209  virtual void printIteration(int ioflag, doublereal damp, size_t label_d, size_t label_t,
210  doublereal inv_t, doublereal t_real, int iter,
211  doublereal update_norm, doublereal resid_norm,
212  doublereal netProdRate[], doublereal CSolnSP[],
213  doublereal resid[],
214  doublereal wtSpecies[], size_t dim, bool do_time);
215
216  //! Print a summary of the solution
217  virtual void printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t,
218  doublereal inv_t, doublereal t_real, int iter,
219  doublereal update_norm, doublereal resid_norm,
220  doublereal netProdRateKinSpecies[], const doublereal CSolnSP[],
221  const doublereal resid[],
222  const doublereal wtSpecies[], const doublereal wtRes[],
223  size_t dim, bool do_time);
224
225  //! Calculate a conservative delta T to use in a pseudo-steady state
226  //! algorithm
227  /*!
228  * This routine calculates a pretty conservative 1/del_t based
229  * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees
230  * diagonal dominance.
231  *
232  * Small surface fractions are allowed to intervene in the del_t
233  * determination, no matter how small. This may be changed.
234  * Now minimum changed to 1.0e-12,
235  *
236  * Maximum time step set to time_scale.
237  *
238  * @param netProdRateSolnSP Output variable. Net production rate
239  * of all of the species in the solution vector.
240  * @param Csoln output variable.
241  * Mole fraction of all of the species in the solution vector
242  * @param label Output variable. Pointer to the value of the
243  * species index (kindexSP) that is controlling
244  * the time step
245  * @param label_old Output variable. Pointer to the value of the
246  * species index (kindexSP) that controlled
247  * the time step at the previous iteration
248  * @param label_factor Output variable. Pointer to the current
249  * factor that is used to indicate the same species
250  * is controlling the time step.
251  *
252  * @param ioflag Level of the output requested.
253  *
254  * @return Returns the 1. / delta T to be used on the next step
255  */
256  virtual doublereal calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
257  size_t* label, size_t* label_old,
258  doublereal* label_factor, int ioflag);
259
260  //! Calculate the solution and residual weights
261  /*!
262  * Calculate the weighting factors for norms wrt both the species
263  * concentration unknowns and the residual unknowns.
264  * @param wtSpecies Weights to use for the soln unknowns. These
265  * are in concentration units
266  * @param wtResid Weights to sue for the residual unknowns.
267  *
268  * @param CSolnSP Solution vector for the surface problem
269  */
270  virtual void calcWeights(doublereal wtSpecies[], doublereal wtResid[],
271  const doublereal CSolnSP[]);
272
273 #ifdef DEBUG_SOLVEPROB
274  //! Utility routine to print a header for high lvls of debugging
275  /*!
276  * @param ioflag Lvl of debugging
277  * @param damp lvl of damping
278  * @param inv_t Inverse of the value of delta T
279  * @param t_real Value of the time
280  * @param iter Iteration number
281  * @param do_time boolean indicating whether time stepping is taking
282  * place
283  */
284  virtual void printIterationHeader(int ioflag, doublereal damp,
285  doublereal inv_t, doublereal t_real, int iter,
286  bool do_time);
287 #endif
288
289  //! Main Function evaluation
290  /*!
291  * This calculates the net production rates of all species
292  *
293  * @param resid output Vector of residuals, length = m_neq
294  * @param CSolnSP Vector of species concentrations, unknowns in the
295  * problem, length = m_neq
296  * @param CSolnSPOld Old Vector of species concentrations, unknowns in the
297  * problem, length = m_neq
298  * @param do_time Calculate a time dependent residual
299  * @param deltaT Delta time for time dependent problem.
300  */
301  virtual void fun_eval(doublereal* const resid, const doublereal* const CSolnSP,
302  const doublereal* const CSolnSPOld, const bool do_time, const doublereal deltaT);
303
304  //! Main routine that calculates the current residual and Jacobian
305  /*!
306  * @param JacCol Vector of pointers to the tops of columns of the
307  * Jacobian to be evaluated.
308  * @param resid output Vector of residuals, length = m_neq
309  * @param CSolnSP Vector of species concentrations, unknowns in the
310  * problem, length = m_neq. These are tweaked in order
311  * to derive the columns of the jacobian.
312  * @param CSolnSPOld Old Vector of species concentrations, unknowns in the
313  * problem, length = m_neq
314  * @param do_time Calculate a time dependent residual
315  * @param deltaT Delta time for time dependent problem.
316  */
317  virtual void resjac_eval(std::vector<doublereal*>& JacCol, doublereal* resid,
318  doublereal* CSolnSP,
319  const doublereal* CSolnSPOld, const bool do_time,
320  const doublereal deltaT);
321
322  //! This function calculates a damping factor for the Newton iteration update
323  //! vector, dxneg, to insure that all solution components stay within prescribed bounds
324  /*!
325  * The default for this class is that all solution components are bounded between zero and one.
326  * this is because the original unknowns were mole fractions and surface site fractions.
327  *
328  * dxneg[] = negative of the update vector.
329  *
330  * The constant "APPROACH" sets the fraction of the distance to the boundary
331  * that the step can take. If the full step would not force any fraction
332  * outside of the bounds, then Newton's method is mostly allowed to operate normally.
333  * There is also some solution damping employed.
334  *
335  * @param x Vector of the current solution components
336  * @param dxneg Vector of the negative of the full solution update vector.
337  * @param dim Size of the solution vector
338  * @param label return int, stating which solution component caused the most damping.
339  */
340  virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label);
341
342  //! residual function pointer to be solved.
344
345  //! Total number of equations to solve in the implicit problem.
346  /*!
347  * Note, this can be zero, and frequently is
348  */
349  size_t m_neq;
350
351  //! m_atol is the absolute tolerance in real units.
353
354  //! m_rtol is the relative error tolerance.
355  doublereal m_rtol;
356
357  //! maximum value of the time step
358  /*!
359  * units = seconds
360  */
361  doublereal m_maxstep;
362
363  //! Temporary vector with length MAX(1, m_neq)
365
366  //! Temporary vector with length MAX(1, m_neq)
368
369  //! Temporary vector with length MAX(1, m_neq)
371
372  //! Temporary vector with length MAX(1, m_neq)
374
375  //! Solution vector
376  /*!
377  * length MAX(1, m_neq)
378  */
380
381  //! Saved initial solution vector
382  /*!
383  * length MAX(1, m_neq)
384  */
386
387  //! Saved solution vector at the old time step
388  /*!
389  * length MAX(1, m_neq)
390  */
392
393  //! Weights for the residual norm calculation
394  /*!
395  * length MAX(1, m_neq)
396  */
398
399  //! Weights for the species concentrations norm calculation
400  /*!
401  * length MAX(1, m_neq)
402  */
404
405  //! Residual for the surface problem
406  /*!
407  * The residual vector of length "dim" that, that has the value
408  * of "sdot" for surface species. The residuals for the bulk
409  * species are a function of the sdots for all species in the bulk
410  * phase. The last residual of each phase enforces {Sum(fractions)
411  * = 1}. After linear solve (dgetrf_ & dgetrs_), resid holds the
412  * update vector.
413  *
414  * length MAX(1, m_neq)
415  */
417
418  //! pivots
419  /*!
420  * length MAX(1, m_neq)
421  */
423
424  //! Vector of pointers to the top of the columns of the jacobians
425  /*!
426  * The "dim" by "dim" computed Jacobian matrix for the
427  * local Newton's method.
428  */
429  std::vector<doublereal*> m_JacCol;
430
431  //! Jacobian
432  /*!
433  * m_neq by m_neq computed Jacobian matrix for the local Newton's method.
434  */
436
437  //! Top bounds for the solution vector
438  /*!
439  * This defaults to 1.0
440  */
442
443  //! Bottom bounds for the solution vector
444  /*!
445  * This defaults to 0.0
446  */
448
449 public:
450  int m_ioflag;
451 };
452 }
453 #endif
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[])
Set the bottom and top bounds on the solution vector.
Definition: solveProb.cpp:609
const int SOLVEPROB_INITIALIZE
Solution Methods.
Definition: solveProb.h:52
virtual doublereal calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[], size_t *label, size_t *label_old, doublereal *label_factor, int ioflag)
Calculate a conservative delta T to use in a pseudo-steady state algorithm.
Definition: solveProb.cpp:557
virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t *label)
This function calculates a damping factor for the Newton iteration update vector, dxneg...
Definition: solveProb.cpp:440
solveProb(ResidEval *resid)
Constructor for the object.
Definition: solveProb.cpp:27
vector_fp m_wtSpecies
Weights for the species concentrations norm calculation.
Definition: solveProb.h:403
int solve(int ifunc, doublereal time_scale, doublereal reltol)
Main routine that actually calculates the pseudo steady state of the surface problem.
Definition: solveProb.cpp:67
virtual void printIteration(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRate[], doublereal CSolnSP[], doublereal resid[], doublereal wtSpecies[], size_t dim, bool do_time)
Printing routine that gets called after every iteration.
Definition: solveProb.cpp:723
solveProb & operator=(const solveProb &right)
Unimplemented private assignment operator.
size_t m_neq
Total number of equations to solve in the implicit problem.
Definition: solveProb.h:349
vector_fp m_numEqn1
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:367
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:33
ResidEval * m_residFunc
residual function pointer to be solved.
Definition: solveProb.h:343
vector_fp m_botBounds
Bottom bounds for the solution vector.
Definition: solveProb.h:447
vector_fp m_CSolnSPOld
Saved solution vector at the old time step.
Definition: solveProb.h:391
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:167
virtual void reportState(doublereal *const CSoln) const
Report the current state of the solution.
Definition: solveProb.cpp:386
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:33
doublereal m_rtol
m_rtol is the relative error tolerance.
Definition: solveProb.h:355
vector_fp m_atol
m_atol is the absolute tolerance in real units.
Definition: solveProb.h:352
doublereal m_maxstep
maximum value of the time step
Definition: solveProb.h:361
vector_fp m_resid
Residual for the surface problem.
Definition: solveProb.h:416
vector_fp m_CSolnSP
Solution vector.
Definition: solveProb.h:379
Array2D m_Jac
Jacobian.
Definition: solveProb.h:435
virtual void fun_eval(doublereal *const resid, const doublereal *const CSolnSP, const doublereal *const CSolnSPOld, const bool do_time, const doublereal deltaT)
Main Function evaluation.
Definition: solveProb.cpp:391
virtual void resjac_eval(std::vector< doublereal * > &JacCol, doublereal *resid, doublereal *CSolnSP, const doublereal *CSolnSPOld, const bool do_time, const doublereal deltaT)
Main routine that calculates the current residual and Jacobian.
Definition: solveProb.cpp:406
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:165
std::vector< doublereal * > m_JacCol
Vector of pointers to the top of the columns of the jacobians.
Definition: solveProb.h:429
vector_fp m_CSolnSPInit
Saved initial solution vector.
Definition: solveProb.h:385
virtual void calcWeights(doublereal wtSpecies[], doublereal wtResid[], const doublereal CSolnSP[])
Calculate the solution and residual weights.
Definition: solveProb.cpp:532
vector_fp m_netProductionRatesSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:364
vector_int m_ipiv
pivots
Definition: solveProb.h:422
virtual void print_header(int ioflag, int ifunc, doublereal time_scale, doublereal reltol, doublereal netProdRate[])
Printing routine that gets called at the start of every invocation.
Definition: solveProb.cpp:626
vector_fp m_CSolnSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:373
vector_fp m_numEqn2
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:370
virtual void printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRateKinSpecies[], const doublereal CSolnSP[], const doublereal resid[], const doublereal wtSpecies[], const doublereal wtRes[], size_t dim, bool do_time)
Print a summary of the solution.
Definition: solveProb.cpp:812
vector_fp m_wtResid
Weights for the residual norm calculation.
Definition: solveProb.h:397
vector_fp m_topBounds
Top bounds for the solution vector.
Definition: solveProb.h:441
Method to solve a pseudo steady state of a nonlinear problem.
Definition: solveProb.h:138