Cantera  3.1.0a1
Loading...
Searching...
No Matches
solveSP.h
Go to the documentation of this file.
1/**
2 * @file solveSP.h Header file for implicit surface problem solver (see @ref
3 * chemkinetics and class @link Cantera::solveSP solveSP@endlink).
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
9#ifndef SOLVESP_H
10#define SOLVESP_H
11
14
15//! @defgroup solvesp_methods Surface Problem Solver Methods
16//! @ingroup surfSolverGroup
17//! @{
18
19//! This assumes that the initial guess supplied to the routine is far from
20//! the correct one. Substantial work plus transient time-stepping is to be
21//! expected to find a solution.
22const int SFLUX_INITIALIZE = 1;
23
24//! Need to solve the surface problem in order to calculate the surface fluxes
25//! of gas-phase species. (Can expect a moderate change in the solution
26//! vector; try to solve the system by direct methods with no damping first,
27//! then try time-stepping if the first method fails). A "time_scale" supplied
28//! here is used in the algorithm to determine when to shut off time-stepping.
29const int SFLUX_RESIDUAL = 2;
30
31//! Calculation of the surface problem is due to the need for a numerical
32//! Jacobian for the gas-problem. The solution is expected to be very close to
33//! the initial guess, and accuracy is needed because solution variables have
34//! been perturbed from nominal values to create Jacobian entries.
35const int SFLUX_JACOBIAN = 3;
36
37//! The transient calculation is performed here for an amount of time
38//! specified by "time_scale". It is not guaranteed to be time-accurate -
39//! just stable and fairly fast. The solution after del_t time is returned,
40//! whether it's converged to a steady state or not. This is a poor man's time
41//! stepping algorithm.
42const int SFLUX_TRANSIENT = 4;
43//! @} End of solvesp_method
44
45//! @defgroup solvesp_bulkFunc Surface Problem Bulk Phase Mode
46//! Functionality expected from the bulk phase. This changes the equations
47//! that will be used to solve for the bulk mole fractions.
48//! @ingroup surfSolverGroup
49//! @{
50
51//! Deposition of a bulk phase is to be expected. Bulk mole fractions are
52//! determined from ratios of growth rates of bulk species.
53const int BULK_DEPOSITION = 1;
54
55//! Etching of a bulk phase is to be expected. Bulk mole fractions are assumed
56//! constant, and given by the initial conditions. This is also used whenever
57//! the condensed phase is part of the larger solution.
58const int BULK_ETCH = 2;
59//! @} End of solvesp_bulkFunc
60
61namespace Cantera
62{
63
64//! @addtogroup surfSolverGroup
65//! @{
66
67//! Method to solve a pseudo steady state surface problem
68/*!
69 * The following class handles solving the surface problem. The calculation
70 * uses Newton's method to obtain the surface fractions of the surface and
71 * bulk species by requiring that the surface species production rate = 0 and
72 * that the either the bulk fractions are proportional to their production
73 * rates or they are constants.
74 *
75 * Currently, the bulk mole fractions are treated as constants. Implementation
76 * of their being added to the unknown solution vector is delayed.
77 *
78 * Lets introduce the unknown vector for the "surface problem". The surface
79 * problem is defined as the evaluation of the surface site fractions for
80 * multiple surface phases. The unknown vector will consist of the vector of
81 * surface concentrations for each species in each surface vector. Species
82 * are grouped first by their surface phases
83 *
84 * - C_i_j = Concentration of the ith species in the jth surface phase
85 * - Nj = number of surface species in the jth surface phase
86 *
87 * The unknown solution vector is defined as follows:
88 *
89 * C_i_j | kindexSP
90 * --------- | ----------
91 * C_0_0 | 0
92 * C_1_0 | 1
93 * C_2_0 | 2
94 * . . . | ...
95 * C_N0-1_0 | N0-1
96 * C_0_1 | N0
97 * C_1_1 | N0+1
98 * C_2_1 | N0+2
99 * . . . | ...
100 * C_N1-1_1 | NO+N1-1
101 *
102 * Note there are a couple of different types of species indices floating
103 * around in the formulation of this object.
104 *
105 * kindexSP: This is the species index in the contiguous vector of unknowns
106 * for the surface problem.
107 *
108 * Note, in the future, BULK_DEPOSITION systems will be added, and the
109 * solveSP unknown vector will get more complicated. It will include the mole
110 * fraction and growth rates of specified bulk phases
111 *
112 * Indices which relate to individual kinetics objects use the suffix KSI
113 * (kinetics species index).
114 *
115 * ## Solution Method
116 *
117 * This routine is typically used within a residual calculation in a large code.
118 * It's typically invoked millions of times for large calculations, and it must
119 * work every time. Therefore, requirements demand that it be robust but also
120 * efficient.
121 *
122 * The solution methodology is largely determined by the `ifunc` parameter,
123 * that is input to the solution object. This parameter may have one of the
124 * values defined in @ref solvesp_methods.
125 *
126 * ### Pseudo time stepping algorithm:
127 * The time step is determined from sdot[], so that the time step
128 * doesn't ever change the value of a variable by more than 100%.
129 *
130 * This algorithm does use a damped Newton's method to relax the equations.
131 * Damping is based on a "delta damping" technique. The solution unknowns
132 * are not allowed to vary too much between iterations.
133 *
134 * `EXTRA_ACCURACY`: A constant that is the ratio of the required update norm
135 * in this Newton iteration compared to that in the nonlinear solver. A value
136 * of 0.1 is used so surface species are safely overconverged.
137 *
138 * Functions called:
139 * - `ct_dgetrf` -- First half of LAPACK direct solve of a full Matrix
140 * - `ct_dgetrs` -- Second half of LAPACK direct solve of a full matrix.
141 * Returns solution vector in the right-hand-side vector, resid.
142 */
144{
145public:
146 //! Constructor for the object
147 /*!
148 * @param surfChemPtr Pointer to the ImplicitSurfChem object that
149 * defines the surface problem to be solved.
150 * @param bulkFunc Integer representing how the bulk phases should be
151 * handled. See @ref solvesp_bulkFunc. Currently,
152 * only the default value of BULK_ETCH is supported.
153 */
154 solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc = BULK_ETCH);
155
156 //! Destructor
157 ~solveSP() = default;
158
159 solveSP(const solveSP&) = delete;
160 solveSP& operator=(const solveSP&) = delete;
161
162public:
163 //! Main routine that actually calculates the pseudo steady state
164 //! of the surface problem
165 /*!
166 * The actual converged solution is returned as part of the internal state
167 * of the InterfaceKinetics objects.
168 *
169 * Uses Newton's method to get the surface fractions of the surface and
170 * bulk species by requiring that the surface species production rate = 0
171 * and that the bulk fractions are proportional to their production rates.
172 *
173 * @param ifunc Determines the type of solution algorithm to be used. See
174 * @ref solvesp_methods for possible values.
175 * @param time_scale Time over which to integrate the surface equations,
176 * where applicable
177 * @param TKelvin Temperature (kelvin)
178 * @param PGas Pressure (pascals)
179 * @param reltol Relative tolerance to use
180 * @param abstol absolute tolerance.
181 * @return 1 if the surface problem is successfully solved.
182 * -1 if the surface problem wasn't solved successfully.
183 * Note the actual converged solution is returned as part of the
184 * internal state of the InterfaceKinetics objects.
185 */
186 int solveSurfProb(int ifunc, double time_scale, double TKelvin,
187 double PGas, double reltol, double abstol);
188
189private:
190 //! Printing routine that optionally gets called at the start of every
191 //! invocation
192 void print_header(int ioflag, int ifunc, double time_scale,
193 int damping, double reltol, double abstol);
194
195 //! Printing routine that gets called after every iteration
196 void printIteration(int ioflag, double damp, int label_d, int label_t,
197 double inv_t, double t_real, size_t iter,
198 double update_norm, double resid_norm,
199 bool do_time, bool final=false);
200
201 //! Calculate a conservative delta T to use in a pseudo-steady state
202 //! algorithm
203 /*!
204 * This routine calculates a pretty conservative 1/del_t based on
205 * MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees diagonal dominance.
206 *
207 * Small surface fractions are allowed to intervene in the del_t
208 * determination, no matter how small. This may be changed.
209 * Now minimum changed to 1.0e-12,
210 *
211 * Maximum time step set to time_scale.
212 *
213 * @param netProdRateSolnSP Output variable. Net production rate of all
214 * of the species in the solution vector.
215 * @param XMolSolnSP output variable. Mole fraction of all of the species
216 * in the solution vector
217 * @param label Output variable. Pointer to the value of the species
218 * index (kindexSP) that is controlling the time step
219 * @param label_old Output variable. Pointer to the value of the species
220 * index (kindexSP) that controlled the time step at the previous
221 * iteration
222 * @param label_factor Output variable. Pointer to the current factor
223 * that is used to indicate the same species is controlling the time
224 * step.
225 * @param ioflag Level of the output requested.
226 * @returns the 1. / delta T to be used on the next step
227 */
228 double calc_t(double netProdRateSolnSP[], double XMolSolnSP[],
229 int* label, int* label_old,
230 double* label_factor, int ioflag);
231
232 //! Calculate the solution and residual weights
233 /*!
234 * @param wtSpecies Weights to use for the soln unknowns. These are in
235 * concentration units
236 * @param wtResid Weights to sue for the residual unknowns.
237 * @param Jac Jacobian. Row sum scaling is used for the Jacobian
238 * @param CSolnSP Solution vector for the surface problem
239 * @param abstol Absolute error tolerance
240 * @param reltol Relative error tolerance
241 */
242 void calcWeights(double wtSpecies[], double wtResid[],
243 const Array2D& Jac, const double CSolnSP[],
244 const double abstol, const double reltol);
245
246 /**
247 * Update the surface states of the surface phases.
248 */
249 void updateState(const double* cSurfSpec);
250
251 //! Update mole fraction vector consisting of unknowns in surface problem
252 /*!
253 * @param XMolSolnSP Vector of mole fractions for the unknowns in the
254 * surface problem.
255 */
256 void updateMFSolnSP(double* XMolSolnSP);
257
258 //! Update the mole fraction vector for a specific kinetic species vector
259 //! corresponding to one InterfaceKinetics object
260 /*!
261 * @param XMolKinSp Mole fraction vector corresponding to a particular
262 * kinetic species for a single InterfaceKinetics Object
263 * This is a vector over all the species in all of the
264 * phases in the InterfaceKinetics object
265 * @param isp ID of the InterfaceKinetics Object.
266 */
267 void updateMFKinSpecies(double* XMolKinSp, int isp);
268
269 //! Update the vector that keeps track of the largest species in each
270 //! surface phase.
271 /*!
272 * @param CSolnSP Vector of the current values of the surface concentrations
273 * in all of the surface species.
274 */
275 void evalSurfLarge(const double* CSolnSP);
276
277 //! Main Function evaluation
278 /*!
279 * @param resid output Vector of residuals, length = m_neq
280 * @param CSolnSP Vector of species concentrations, unknowns in the
281 * problem, length = m_neq
282 * @param CSolnOldSP Old Vector of species concentrations, unknowns in the
283 * problem, length = m_neq
284 * @param do_time Calculate a time dependent residual
285 * @param deltaT Delta time for time dependent problem.
286 */
287 void fun_eval(double* resid, const double* CSolnSP,
288 const double* CSolnOldSP, const bool do_time, const double deltaT);
289
290 //! Main routine that calculates the current residual and Jacobian
291 /*!
292 * @param jac Jacobian to be evaluated.
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. These are tweaked in order
296 * to derive the columns of the Jacobian.
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 void resjac_eval(DenseMatrix& jac, double* resid,
303 double* CSolnSP,
304 const double* CSolnSPOld, const bool do_time,
305 const double deltaT);
306
307 //! Vector of interface kinetics objects
308 /*!
309 * Each of these is associated with one and only one surface phase.
310 */
311 vector<InterfaceKinetics*> &m_objects;
312
313 //! Total number of equations to solve in the implicit problem.
314 /*!
315 * Note, this can be zero, and frequently is
316 */
317 size_t m_neq = 0;
318
319 //! This variable determines how the bulk phases are to be handled
320 /*!
321 * Possible values are given in @ref solvesp_bulkFunc.
322 */
324
325 //! Number of surface phases in the surface problem
326 /*!
327 * This number is equal to the number of InterfaceKinetics objects
328 * in the problem. (until further noted)
329 */
330 size_t m_numSurfPhases = 0;
331
332 //! Total number of surface species in all surface phases.
333 /*!
334 * This is also the number of equations to solve for m_mode=0 system
335 * It's equal to the sum of the number of species in each of the
336 * m_numSurfPhases.
337 */
339
340 //! Mapping between the surface phases and the InterfaceKinetics objects
341 /*!
342 * Currently this is defined to be a 1-1 mapping (and probably assumed
343 * in some places)
344 * m_surfKinObjID[i] = i
345 */
347
348 //! Vector of length number of surface phases containing
349 //! the number of surface species in each phase
350 /*!
351 * Length is equal to the number of surface phases, m_numSurfPhases
352 */
353 vector<size_t> m_nSpeciesSurfPhase;
354
355 //! Vector of surface phase pointers
356 /*!
357 * This is created during the constructor
358 * Length is equal to the number of surface phases, m_numSurfPhases
359 */
360 vector<SurfPhase*> m_ptrsSurfPhase;
361
362 //! Index of the start of the unknowns for each solution phase
363 /*!
364 * i_eqn = m_eqnIndexStartPhase[isp]
365 *
366 * isp is the phase id in the list of phases solved by the
367 * surface problem.
368 *
369 * i_eqn is the equation number of the first unknown in the
370 * solution vector corresponding to isp'th phase.
371 */
373
374 //! Total number of volumetric condensed phases included in the steady state
375 //! problem handled by this routine.
376 /*!
377 * This is equal to or less than the total number of volumetric phases in
378 * all of the InterfaceKinetics objects. We usually do not include bulk
379 * phases. Bulk phases are only included in the calculation when their
380 * domain isn't included in the underlying continuum model conservation
381 * equation system.
382 *
383 * This is equal to 0, for the time being
384 */
386
387 //! Vector of number of species in the m_numBulkPhases phases.
388 /*!
389 * Length is number of bulk phases
390 */
391 vector<size_t> m_numBulkSpecies;
392
393 //! Total number of species in all bulk phases.
394 /*!
395 * This is also the number of bulk equations to solve when bulk equation
396 * solving is turned on.
397 */
399
400 //! Vector of bulk phase pointers, length is equal to m_numBulkPhases.
401 vector<ThermoPhase*> m_bulkPhasePtrs;
402
403 //! Index between the equation index and the position in the kinetic
404 //! species array for the appropriate kinetics operator
405 /*!
406 * Length = m_neq.
407 *
408 * ksp = m_kinSpecIndex[ieq]
409 * ksp is the kinetic species index for the ieq'th equation.
410 */
411 vector<size_t> m_kinSpecIndex;
412
413 //! Index between the equation index and the index of the
414 //! InterfaceKinetics object
415 /*!
416 * Length m_neq
417 */
418 vector<size_t> m_kinObjIndex;
419
420 //! Vector containing the indices of the largest species
421 //! in each surface phase
422 /*!
423 * `k = m_spSurfLarge[i]` where `k` is the local species index, that is, it
424 * varies from 0 to (num species in phase - 1) and `i` is the surface
425 * phase index in the problem. Length is equal to #m_numSurfPhases.
426 */
427 vector<size_t> m_spSurfLarge;
428
429 //! Maximum number of species in any single kinetics operator
430 //! -> also maxed wrt the total # of solution species
431 size_t m_maxTotSpecies = 0;
432
433 //! Temporary vector with length equal to max m_maxTotSpecies
435
436 //! Temporary vector with length equal to max m_maxTotSpecies
437 vector<double> m_numEqn1;
438
439 //! Temporary vector with length equal to max m_maxTotSpecies
440 vector<double> m_numEqn2;
441
442 //! Temporary vector with length equal to max m_maxTotSpecies
443 vector<double> m_CSolnSave;
444
445 //! Solution vector. length MAX(1, m_neq)
446 vector<double> m_CSolnSP;
447
448 //! Saved initial solution vector. length MAX(1, m_neq)
449 vector<double> m_CSolnSPInit;
450
451 //! Saved solution vector at the old time step. length MAX(1, m_neq)
452 vector<double> m_CSolnSPOld;
453
454 //! Weights for the residual norm calculation. length MAX(1, m_neq)
455 vector<double> m_wtResid;
456
457 //! Weights for the species concentrations norm calculation
458 /*!
459 * length MAX(1, m_neq)
460 */
461 vector<double> m_wtSpecies;
462
463 //! Residual for the surface problem
464 /*!
465 * The residual vector of length "dim" that, that has the value of "sdot"
466 * for surface species. The residuals for the bulk species are a function
467 * of the sdots for all species in the bulk phase. The last residual of
468 * each phase enforces {Sum(fractions) = 1}. After linear solve (dgetrf_ &
469 * dgetrs_), resid holds the update vector.
470 *
471 * length MAX(1, m_neq)
472 */
473 vector<double> m_resid;
474
475 //! Vector of mole fractions. length m_maxTotSpecies
476 vector<double> m_XMolKinSpecies;
477
478 //! Jacobian. m_neq by m_neq computed Jacobian matrix for the local
479 //! Newton's method.
481
482public:
483 int m_ioflag = 0;
484};
485
486//! @} End of surfSolverGroup
487
488}
489#endif
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
Advances the surface coverages of the associated set of SurfacePhase objects in time.
Method to solve a pseudo steady state surface problem.
Definition solveSP.h:144
vector< size_t > m_numBulkSpecies
Vector of number of species in the m_numBulkPhases phases.
Definition solveSP.h:391
vector< double > m_wtSpecies
Weights for the species concentrations norm calculation.
Definition solveSP.h:461
size_t m_numTotSurfSpecies
Total number of surface species in all surface phases.
Definition solveSP.h:338
void evalSurfLarge(const double *CSolnSP)
Update the vector that keeps track of the largest species in each surface phase.
Definition solveSP.cpp:303
void updateState(const double *cSurfSpec)
Update the surface states of the surface phases.
Definition solveSP.cpp:272
vector< size_t > m_indexKinObjSurfPhase
Mapping between the surface phases and the InterfaceKinetics objects.
Definition solveSP.h:346
vector< size_t > m_kinObjIndex
Index between the equation index and the index of the InterfaceKinetics object.
Definition solveSP.h:418
vector< double > m_resid
Residual for the surface problem.
Definition solveSP.h:473
vector< ThermoPhase * > m_bulkPhasePtrs
Vector of bulk phase pointers, length is equal to m_numBulkPhases.
Definition solveSP.h:401
vector< size_t > m_spSurfLarge
Vector containing the indices of the largest species in each surface phase.
Definition solveSP.h:427
vector< double > m_netProductionRatesSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition solveSP.h:434
double calc_t(double netProdRateSolnSP[], double XMolSolnSP[], int *label, int *label_old, double *label_factor, int ioflag)
Calculate a conservative delta T to use in a pseudo-steady state algorithm.
Definition solveSP.cpp:563
~solveSP()=default
Destructor.
int m_bulkFunc
This variable determines how the bulk phases are to be handled.
Definition solveSP.h:323
vector< double > m_wtResid
Weights for the residual norm calculation. length MAX(1, m_neq)
Definition solveSP.h:455
vector< InterfaceKinetics * > & m_objects
Vector of interface kinetics objects.
Definition solveSP.h:311
void fun_eval(double *resid, const double *CSolnSP, const double *CSolnOldSP, const bool do_time, const double deltaT)
Main Function evaluation.
Definition solveSP.cpp:319
vector< double > m_CSolnSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition solveSP.h:443
void updateMFSolnSP(double *XMolSolnSP)
Update mole fraction vector consisting of unknowns in surface problem.
Definition solveSP.cpp:286
vector< double > m_CSolnSPOld
Saved solution vector at the old time step. length MAX(1, m_neq)
Definition solveSP.h:452
void printIteration(int ioflag, double damp, int label_d, int label_t, double inv_t, double t_real, size_t iter, double update_norm, double resid_norm, bool do_time, bool final=false)
Printing routine that gets called after every iteration.
Definition solveSP.cpp:651
vector< double > m_CSolnSP
Solution vector. length MAX(1, m_neq)
Definition solveSP.h:446
vector< size_t > m_kinSpecIndex
Index between the equation index and the position in the kinetic species array for the appropriate ki...
Definition solveSP.h:411
size_t m_numTotBulkSpeciesSS
Total number of species in all bulk phases.
Definition solveSP.h:398
vector< double > m_numEqn1
Temporary vector with length equal to max m_maxTotSpecies.
Definition solveSP.h:437
vector< double > m_CSolnSPInit
Saved initial solution vector. length MAX(1, m_neq)
Definition solveSP.h:449
vector< double > m_XMolKinSpecies
Vector of mole fractions. length m_maxTotSpecies.
Definition solveSP.h:476
vector< double > m_numEqn2
Temporary vector with length equal to max m_maxTotSpecies.
Definition solveSP.h:440
vector< size_t > m_eqnIndexStartSolnPhase
Index of the start of the unknowns for each solution phase.
Definition solveSP.h:372
vector< size_t > m_nSpeciesSurfPhase
Vector of length number of surface phases containing the number of surface species in each phase.
Definition solveSP.h:353
size_t m_numSurfPhases
Number of surface phases in the surface problem.
Definition solveSP.h:330
void print_header(int ioflag, int ifunc, double time_scale, int damping, double reltol, double abstol)
Printing routine that optionally gets called at the start of every invocation.
Definition solveSP.cpp:602
void updateMFKinSpecies(double *XMolKinSp, int isp)
Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceK...
Definition solveSP.cpp:294
DenseMatrix m_Jac
Jacobian.
Definition solveSP.h:480
size_t m_numBulkPhasesSS
Total number of volumetric condensed phases included in the steady state problem handled by this rout...
Definition solveSP.h:385
size_t m_neq
Total number of equations to solve in the implicit problem.
Definition solveSP.h:317
vector< SurfPhase * > m_ptrsSurfPhase
Vector of surface phase pointers.
Definition solveSP.h:360
size_t m_maxTotSpecies
Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution s...
Definition solveSP.h:431
int solveSurfProb(int ifunc, double time_scale, double TKelvin, double PGas, double reltol, double abstol)
Main routine that actually calculates the pseudo steady state of the surface problem.
Definition solveSP.cpp:89
void calcWeights(double wtSpecies[], double wtResid[], const Array2D &Jac, const double CSolnSP[], const double abstol, const double reltol)
Calculate the solution and residual weights.
Definition solveSP.cpp:530
void resjac_eval(DenseMatrix &jac, double *resid, double *CSolnSP, const double *CSolnSPOld, const bool do_time, const double deltaT)
Main routine that calculates the current residual and Jacobian.
Definition solveSP.cpp:418
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition solveSP.h:58
const int BULK_DEPOSITION
Deposition of a bulk phase is to be expected.
Definition solveSP.h:53
const int SFLUX_TRANSIENT
The transient calculation is performed here for an amount of time specified by "time_scale".
Definition solveSP.h:42
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species.
Definition solveSP.h:29
const int SFLUX_JACOBIAN
Calculation of the surface problem is due to the need for a numerical Jacobian for the gas-problem.
Definition solveSP.h:35
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one.
Definition solveSP.h:22
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564