Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
4  * Cantera::solveProb solveProb\endlink).
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 
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  * @deprecated Unused. To be removed after Cantera 2.2.
138  */
140 {
141 public:
142 
143  //! Constructor for the object
144  solveProb(ResidEval* resid);
145 
146  virtual ~solveProb() {}
147 
148 private:
149 
150  //! Unimplemented private copy constructor
151  solveProb(const solveProb& right);
152 
153  //! Unimplemented private assignment operator
154  solveProb& operator=(const solveProb& right);
155 
156 public:
157  //! Main routine that actually calculates the pseudo steady state
158  //! of the surface problem
159  /*!
160  * The actual converged solution is returned as part of the
161  * internal state of the InterfaceKinetics objects.
162  *
163  * @param ifunc Determines the type of solution algorithm to be
164  * used. Possible values are SOLVEPROB_INITIALIZE ,
165  * SOLVEPROB_RESIDUAL SOLVEPROB_JACOBIAN SOLVEPROB_TRANSIENT .
166  *
167  * @param time_scale Time over which to integrate the surface equations,
168  * where applicable
169  *
170  * @param reltol Relative tolerance to use
171  *
172  * @return Returns 1 if the surface problem is successfully solved.
173  * Returns -1 if the surface problem wasn't solved successfully.
174  * Note the actual converged solution is returned as part of the
175  * internal state of the InterfaceKinetics objects.
176  */
177  int solve(int ifunc, doublereal time_scale, doublereal reltol);
178 
179  //! Report the current state of the solution
180  /*!
181  * @param[out] CSoln solution vector for the nonlinear problem
182  */
183  virtual void reportState(doublereal* const CSoln) const;
184 
185  //! Set the bottom and top bounds on the solution vector
186  /*!
187  * The default is for the bottom is 0.0, while the default for the top is 1.0
188  *
189  * @param botBounds Vector of bottom bounds
190  * @param topBounds vector of top bounds
191  */
192  virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[]);
193 
194  void setAtol(const doublereal atol[]);
195  void setAtolConst(const doublereal atolconst);
196 
197 private:
198  //! Printing routine that gets called at the start of every invocation
199  virtual void print_header(int ioflag, int ifunc, doublereal time_scale,
200  doublereal reltol,
201  doublereal netProdRate[]);
202 
203 #ifdef DEBUG_SOLVEPROB
204  //! Prints out the residual and Jacobian
205  virtual void printResJac(int ioflag, int neq, const Array2D& Jac,
206  doublereal resid[], doublereal wtResid[], doublereal norm);
207 #endif
208 
209  //! Printing routine that gets called after every iteration
210  virtual void printIteration(int ioflag, doublereal damp, size_t label_d, size_t label_t,
211  doublereal inv_t, doublereal t_real, int iter,
212  doublereal update_norm, doublereal resid_norm,
213  doublereal netProdRate[], doublereal CSolnSP[],
214  doublereal resid[],
215  doublereal wtSpecies[], size_t dim, bool do_time);
216 
217  //! Print a summary of the solution
218  virtual void printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t,
219  doublereal inv_t, doublereal t_real, int iter,
220  doublereal update_norm, doublereal resid_norm,
221  doublereal netProdRateKinSpecies[], const doublereal CSolnSP[],
222  const doublereal resid[],
223  const doublereal wtSpecies[], const doublereal wtRes[],
224  size_t dim, bool do_time);
225 
226  //! Calculate a conservative delta T to use in a pseudo-steady state
227  //! algorithm
228  /*!
229  * This routine calculates a pretty conservative 1/del_t based
230  * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees
231  * diagonal dominance.
232  *
233  * Small surface fractions are allowed to intervene in the del_t
234  * determination, no matter how small. This may be changed.
235  * Now minimum changed to 1.0e-12,
236  *
237  * Maximum time step set to time_scale.
238  *
239  * @param netProdRateSolnSP Output variable. Net production rate
240  * of all of the species in the solution vector.
241  * @param Csoln output variable.
242  * Mole fraction of all of the species in the solution vector
243  * @param label Output variable. Pointer to the value of the
244  * species index (kindexSP) that is controlling
245  * the time step
246  * @param label_old Output variable. Pointer to the value of the
247  * species index (kindexSP) that controlled
248  * the time step at the previous iteration
249  * @param label_factor Output variable. Pointer to the current
250  * factor that is used to indicate the same species
251  * is controlling the time step.
252  *
253  * @param ioflag Level of the output requested.
254  *
255  * @return Returns the 1. / delta T to be used on the next step
256  */
257  virtual doublereal calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
258  size_t* label, size_t* label_old,
259  doublereal* label_factor, int ioflag);
260 
261  //! Calculate the solution and residual weights
262  /*!
263  * Calculate the weighting factors for norms wrt both the species
264  * concentration unknowns and the residual unknowns.
265  * @param wtSpecies Weights to use for the soln unknowns. These
266  * are in concentration units
267  * @param wtResid Weights to sue for the residual unknowns.
268  *
269  * @param CSolnSP Solution vector for the surface problem
270  */
271  virtual void calcWeights(doublereal wtSpecies[], doublereal wtResid[],
272  const doublereal CSolnSP[]);
273 
274 #ifdef DEBUG_SOLVEPROB
275  //! Utility routine to print a header for high lvls of debugging
276  /*!
277  * @param ioflag Lvl of debugging
278  * @param damp lvl of damping
279  * @param inv_t Inverse of the value of delta T
280  * @param t_real Value of the time
281  * @param iter Iteration number
282  * @param do_time boolean indicating whether time stepping is taking
283  * place
284  */
285  virtual void printIterationHeader(int ioflag, doublereal damp,
286  doublereal inv_t, doublereal t_real, int iter,
287  bool do_time);
288 #endif
289 
290  //! Main Function evaluation
291  /*!
292  * This calculates the net production rates of all species
293  *
294  * @param resid output Vector of residuals, length = m_neq
295  * @param CSolnSP Vector of species concentrations, unknowns in the
296  * problem, length = m_neq
297  * @param CSolnSPOld Old Vector of species concentrations, unknowns in the
298  * problem, length = m_neq
299  * @param do_time Calculate a time dependent residual
300  * @param deltaT Delta time for time dependent problem.
301  */
302  virtual void fun_eval(doublereal* const resid, const doublereal* const CSolnSP,
303  const doublereal* const CSolnSPOld, const bool do_time, const doublereal deltaT);
304 
305  //! Main routine that calculates the current residual and Jacobian
306  /*!
307  * @param JacCol Vector of pointers to the tops of columns of the
308  * Jacobian to be evaluated.
309  * @param resid output Vector of residuals, length = m_neq
310  * @param CSolnSP Vector of species concentrations, unknowns in the
311  * problem, length = m_neq. These are tweaked in order
312  * to derive the columns of the Jacobian.
313  * @param CSolnSPOld Old Vector of species concentrations, unknowns in the
314  * problem, length = m_neq
315  * @param do_time Calculate a time dependent residual
316  * @param deltaT Delta time for time dependent problem.
317  */
318  virtual void resjac_eval(std::vector<doublereal*>& JacCol, doublereal* resid,
319  doublereal* CSolnSP,
320  const doublereal* CSolnSPOld, const bool do_time,
321  const doublereal deltaT);
322 
323  //! This function calculates a damping factor for the Newton iteration update
324  //! vector, dxneg, to insure that all solution components stay within prescribed bounds
325  /*!
326  * The default for this class is that all solution components are bounded between zero and one.
327  * this is because the original unknowns were mole fractions and surface site fractions.
328  *
329  * dxneg[] = negative of the update vector.
330  *
331  * The constant "APPROACH" sets the fraction of the distance to the boundary
332  * that the step can take. If the full step would not force any fraction
333  * outside of the bounds, then Newton's method is mostly allowed to operate normally.
334  * There is also some solution damping employed.
335  *
336  * @param x Vector of the current solution components
337  * @param dxneg Vector of the negative of the full solution update vector.
338  * @param dim Size of the solution vector
339  * @param label return int, stating which solution component caused the most damping.
340  */
341  virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label);
342 
343  //! residual function pointer to be solved.
345 
346  //! Total number of equations to solve in the implicit problem.
347  /*!
348  * Note, this can be zero, and frequently is
349  */
350  size_t m_neq;
351 
352  //! m_atol is the absolute tolerance in real units.
354 
355  //! m_rtol is the relative error tolerance.
356  doublereal m_rtol;
357 
358  //! maximum value of the time step
359  /*!
360  * units = seconds
361  */
362  doublereal m_maxstep;
363 
364  //! Temporary vector with length MAX(1, m_neq)
366 
367  //! Temporary vector with length MAX(1, m_neq)
369 
370  //! Temporary vector with length MAX(1, m_neq)
372 
373  //! Temporary vector with length MAX(1, m_neq)
375 
376  //! Solution vector
377  /*!
378  * length MAX(1, m_neq)
379  */
381 
382  //! Saved initial solution vector
383  /*!
384  * length MAX(1, m_neq)
385  */
387 
388  //! Saved solution vector at the old time step
389  /*!
390  * length MAX(1, m_neq)
391  */
393 
394  //! Weights for the residual norm calculation
395  /*!
396  * length MAX(1, m_neq)
397  */
399 
400  //! Weights for the species concentrations norm calculation
401  /*!
402  * length MAX(1, m_neq)
403  */
405 
406  //! Residual for the surface problem
407  /*!
408  * The residual vector of length "dim" that, that has the value
409  * of "sdot" for surface species. The residuals for the bulk
410  * species are a function of the sdots for all species in the bulk
411  * phase. The last residual of each phase enforces {Sum(fractions)
412  * = 1}. After linear solve (dgetrf_ & dgetrs_), resid holds the
413  * update vector.
414  *
415  * length MAX(1, m_neq)
416  */
418 
419  //! Vector of pointers to the top of the columns of the Jacobians
420  /*!
421  * The "dim" by "dim" computed Jacobian matrix for the
422  * local Newton's method.
423  */
424  std::vector<doublereal*> m_JacCol;
425 
426  //! Jacobian
427  /*!
428  * m_neq by m_neq computed Jacobian matrix for the local Newton's method.
429  */
431 
432  //! Top bounds for the solution vector
433  /*!
434  * This defaults to 1.0
435  */
437 
438  //! Bottom bounds for the solution vector
439  /*!
440  * This defaults to 0.0
441  */
443 
444 public:
445  int m_ioflag;
446 };
447 }
448 #endif
Dense, Square (not sparse) matrices.
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[])
Set the bottom and top bounds on the solution vector.
Definition: solveProb.cpp:589
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:536
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:424
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:404
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:63
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:703
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:350
vector_fp m_numEqn1
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:368
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:29
ResidEval * m_residFunc
residual function pointer to be solved.
Definition: solveProb.h:344
vector_fp m_botBounds
Bottom bounds for the solution vector.
Definition: solveProb.h:442
vector_fp m_CSolnSPOld
Saved solution vector at the old time step.
Definition: solveProb.h:392
virtual void reportState(doublereal *const CSoln) const
Report the current state of the solution.
Definition: solveProb.cpp:370
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:34
doublereal m_rtol
m_rtol is the relative error tolerance.
Definition: solveProb.h:356
vector_fp m_atol
m_atol is the absolute tolerance in real units.
Definition: solveProb.h:353
A class for full (non-sparse) matrices with Fortran-compatible data storage.
Definition: SquareMatrix.h:26
doublereal m_maxstep
maximum value of the time step
Definition: solveProb.h:362
vector_fp m_resid
Residual for the surface problem.
Definition: solveProb.h:417
vector_fp m_CSolnSP
Solution vector.
Definition: solveProb.h:380
SquareMatrix m_Jac
Jacobian.
Definition: solveProb.h:430
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:375
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:390
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
std::vector< doublereal * > m_JacCol
Vector of pointers to the top of the columns of the Jacobians.
Definition: solveProb.h:424
vector_fp m_CSolnSPInit
Saved initial solution vector.
Definition: solveProb.h:386
virtual void calcWeights(doublereal wtSpecies[], doublereal wtResid[], const doublereal CSolnSP[])
Calculate the solution and residual weights.
Definition: solveProb.cpp:512
vector_fp m_netProductionRatesSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:365
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:606
vector_fp m_CSolnSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:374
vector_fp m_numEqn2
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:371
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:791
vector_fp m_wtResid
Weights for the residual norm calculation.
Definition: solveProb.h:398
vector_fp m_topBounds
Top bounds for the solution vector.
Definition: solveProb.h:436
Method to solve a pseudo steady state of a nonlinear problem.
Definition: solveProb.h:139