18 static doublereal calc_damping(doublereal* x, doublereal* dx,
size_t dim,
int*);
19 static doublereal calcWeightedNorm(
const doublereal [],
const doublereal dx[],
size_t);
24 m_SurfChemPtr(surfChemPtr),
25 m_objects(surfChemPtr->getObjects()),
29 m_numTotSurfSpecies(0),
31 m_numTotBulkSpeciesSS(0),
39 for (
size_t n = 0; n <
m_objects.size(); n++) {
42 if (surfPhaseIndex !=
npos) {
48 "InterfaceKinetics object has no surface phase");
53 "Inconsistent ThermoPhase object within " 54 "InterfaceKinetics object");
73 size_t tsp =
m_objects[n]->nTotalSpecies();
96 for (
size_t k = 0; k < nsp; k++, kindexSP++) {
103 size_t dim1 = std::max<size_t>(1,
m_neq);
114 doublereal PGas, doublereal reltol, doublereal abstol)
116 doublereal EXTRA_ACCURACY = 0.001;
118 EXTRA_ACCURACY *= 0.001;
123 doublereal label_factor = 1.0;
126 doublereal deltaT = 1.0E-10;
128 doublereal inv_t = 0.0;
129 doublereal t_real = 0.0, update_norm = 1.0E6;
130 bool do_time =
false, not_converged =
true;
131 m_ioflag = std::min(m_ioflag, 1);
155 print_header(m_ioflag, ifunc, time_scale,
true, reltol, abstol);
160 not_converged =
false;
165 while (not_converged && iter < iter_max) {
185 &label_t, &label_t_old, &label_factor, m_ioflag);
188 }
else if (tmp > 2.0*inv_t) {
196 tmp = t_real + 1.0/inv_t;
197 if (tmp > time_scale) {
198 inv_t = 1.0/(time_scale - t_real);
237 for (
size_t irow = 0; irow <
m_neq; irow++) {
240 for (
size_t irow = 0; irow <
m_neq; irow++) {
246 t_real += damp/inv_t;
250 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
251 update_norm, resid_norm, do_time);
255 not_converged = (t_real < time_scale);
258 if (t_real > time_scale ||
259 (resid_norm < 1.0e-7 &&
260 update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
264 not_converged = ((update_norm > EXTRA_ACCURACY) ||
265 (resid_norm > EXTRA_ACCURACY));
272 if (not_converged && m_ioflag) {
273 writelog(
"#$#$#$# Error in solveSP $#$#$#$ \n");
274 writelogf(
"Newton iter on surface species did not converge, " 275 "update_norm = %e \n", update_norm);
285 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
286 update_norm, resid_norm, do_time,
true);
290 if (update_norm > 1.0) {
316 for (
size_t iph = 0; iph < kin->
nPhases(); iph++) {
326 doublereal Clarge = CSolnSP[kindexSP];
330 if (CSolnSP[kindexSP] > Clarge) {
331 Clarge = CSolnSP[kindexSP];
339 const doublereal* CSolnOld,
const bool do_time,
340 const doublereal deltaT)
343 doublereal lenScale = 1.0E-9;
361 size_t kins = kindexSP;
363 for (k = 0; k < nsp; k++, kindexSP++) {
365 (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
371 resid[kspecial] = sd;
372 for (k = 0; k < nsp; k++) {
373 resid[kspecial] -= CSoln[kins + k];
383 size_t kins = kindexSP;
385 for (k = 0; k < nsp; k++, kindexSP++) {
390 resid[kspecial] = sd;
391 for (k = 0; k < nsp; k++) {
392 resid[kspecial] -= CSoln[kins + k];
406 for (k = 0; k < nsp; k++) {
412 for (k = 0; k < nsp; k++) {
413 resid[kindexSP] -= CSoln[kindexSP + k];
416 for (k = 1; k < nsp; k++) {
418 resid[kindexSP + k] = XBlk[k] * grRate
421 resid[kindexSP + k] = XBlk[k] * grRate;
428 for (k = 1; k < nsp; k++) {
429 resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp);
433 for (k = 1; k < nsp; k++) {
434 resid[kindexSP + k] +=
436 (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]);
446 doublereal resid[], doublereal CSoln[],
447 const doublereal CSolnOld[],
const bool do_time,
448 const doublereal deltaT)
450 size_t kColIndex = 0;
452 fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
457 for (
size_t kCol = 0; kCol < nsp; kCol++) {
458 double cSave = CSoln[kColIndex];
459 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
460 CSoln[kColIndex] += dc;
462 for (
size_t i = 0; i <
m_neq; i++) {
463 jac(i, kColIndex) = (
m_numEqn2[i] - resid[i])/dc;
465 CSoln[kColIndex] = cSave;
474 for (
size_t kCol = 0; kCol < nsp; kCol++) {
475 double cSave = CSoln[kColIndex];
476 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
477 CSoln[kColIndex] += dc;
479 for (
size_t i = 0; i <
m_neq; i++) {
480 jac(i, kColIndex) = (
m_numEqn2[i] - resid[i])/dc;
482 CSoln[kColIndex] = cSave;
500 static doublereal calc_damping(doublereal x[], doublereal dxneg[],
size_t dim,
int* label)
502 const doublereal APPROACH = 0.80;
503 doublereal damp = 1.0;
504 static doublereal damp_old = 1.0;
507 for (
size_t i = 0; i < dim; i++) {
509 double xnew = x[i] - damp * dxneg[i];
514 double xtop = 1.0 - 0.1*fabs(1.0-x[i]);
515 double xbot = fabs(x[i]*0.1) - 1.0e-16;
517 damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
519 }
else if (xnew < xbot) {
520 damp = APPROACH * x[i] / dxneg[i];
522 }
else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
523 damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
527 damp = std::max(damp, 1e-2);
531 if (damp > damp_old*3) {
546 static doublereal calcWeightedNorm(
const doublereal wtX[],
const doublereal dx[],
size_t dim)
548 doublereal norm = 0.0;
552 for (
size_t i = 0; i < dim; i++) {
553 norm += pow(dx[i] / wtX[i], 2);
555 return sqrt(norm/dim);
559 const Array2D& Jac,
const doublereal CSoln[],
560 const doublereal abstol,
const doublereal reltol)
568 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
575 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
583 for (
size_t k = 0; k <
m_neq; k++) {
585 for (
size_t jcol = 0; jcol <
m_neq; jcol++) {
586 wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
592 doublereal XMolSolnSP[],
593 int* label,
int* label_old,
594 doublereal* label_factor,
int ioflag)
596 doublereal inv_timeScale = 1.0E-10;
611 size_t kspindex = kstart + k;
612 netProdRateSolnSP[kindexSP] =
m_numEqn1[kspindex];
613 double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10);
615 tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
616 if (netProdRateSolnSP[kindexSP]> 0.0) {
619 if (tmp > inv_timeScale) {
621 *label = int(kindexSP);
628 if (*label == *label_old) {
629 *label_factor *= 1.5;
634 return inv_timeScale / *label_factor;
638 int damping, doublereal reltol, doublereal abstol)
641 writelog(
"\n================================ SOLVESP CALL SETUP " 642 "========================================\n");
644 writelog(
"\n SOLVESP Called with Initialization turned on\n");
645 writelogf(
" Time scale input = %9.3e\n", time_scale);
647 writelog(
"\n SOLVESP Called to calculate steady state residual\n");
648 writelog(
" from a good initial guess\n");
650 writelog(
"\n SOLVESP Called to calculate steady state Jacobian\n");
651 writelog(
" from a good initial guess\n");
653 writelog(
"\n SOLVESP Called to integrate surface in time\n");
654 writelogf(
" for a total of %9.3e sec\n", time_scale);
657 "Unknown ifunc flag = {}", ifunc);
661 writelog(
" The composition of the Bulk Phases will be calculated\n");
663 writelog(
" Bulk Phases have fixed compositions\n");
675 writelogf(
" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
679 writelog(
"\n\n\t Iter Time Del_t Damp DelX " 680 " Resid Name-Time Name-Damp\n");
681 writelog(
"\t -----------------------------------------------" 682 "------------------------------------\n");
687 int label_t, doublereal inv_t, doublereal t_real,
688 size_t iter, doublereal update_norm,
689 doublereal resid_norm,
bool do_time,
bool final)
698 writelogf(
"%9.4e %9.4e ", t_real, 1.0/inv_t);
700 writeline(
' ', 22,
false);
705 writeline(
' ', 11,
false);
707 writelogf(
"%9.4e %9.4e", update_norm, resid_norm);
713 writeline(
' ', 16,
false);
void printIteration(int ioflag, doublereal damp, int label_d, int label_t, doublereal inv_t, doublereal t_real, size_t iter, doublereal update_norm, doublereal resid_norm, bool do_time, bool final=false)
Printing routine that gets called after every iteration.
int m_bulkFunc
This variable determines how the bulk phases are to be handled.
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
std::vector< InterfaceKinetics * > & m_objects
Vector of interface kinetics objects.
size_t m_neq
Total number of equations to solve in the implicit problem.
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
const size_t npos
index returned by functions to indicate "no position"
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
vector_fp m_CSolnSPInit
Saved initial solution vector. length MAX(1, m_neq)
size_t m_numBulkPhasesSS
Total number of volumetric condensed phases included in the steady state problem handled by this rout...
size_t nSpecies() const
Returns the number of species in the phase.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void updateMFSolnSP(doublereal *XMolSolnSP)
Update mole fraction vector consisting of unknowns in surface problem.
vector_fp m_CSolnSP
Solution vector. length MAX(1, m_neq)
const int BULK_ETCH
Etching of a bulk phase is to be expected.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
void evalSurfLarge(const doublereal *CSolnSP)
Update the vector that keeps track of the largest species in each surface phase.
void resjac_eval(DenseMatrix &jac, doublereal *resid, doublereal *CSolnSP, const doublereal *CSolnSPOld, const bool do_time, const doublereal deltaT)
Main routine that calculates the current residual and Jacobian.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
size_t m_numTotSurfSpecies
Total number of surface species in all surface phases.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
void print_header(int ioflag, int ifunc, doublereal time_scale, int damping, doublereal reltol, doublereal abstol)
Printing routine that optionally gets called at the start of every invocation.
int solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, doublereal PGas, doublereal reltol, doublereal abstol)
Main routine that actually calculates the pseudo steady state of the surface problem.
doublereal molarDensity() const
Molar density (kmol/m^3).
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one...
std::vector< size_t > m_eqnIndexStartSolnPhase
Index of the start of the unknowns for each solution phase.
A kinetics manager for heterogeneous reaction mechanisms.
vector_fp m_wtResid
Weights for the residual norm calculation. length MAX(1, m_neq)
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species...
vector_fp m_XMolKinSpecies
Vector of mole fractions. length m_maxTotSpecies.
doublereal calc_t(doublereal netProdRateSolnSP[], doublereal XMolSolnSP[], int *label, int *label_old, doublereal *label_factor, int ioflag)
Calculate a conservative delta T to use in a pseudo-steady state algorithm.
vector_fp m_numEqn1
Temporary vector with length equal to max m_maxTotSpecies.
const int BULK_DEPOSITION
Deposition of a bulk phase is to be expected.
const int SFLUX_JACOBIAN
Calculation of the surface problem is due to the need for a numerical Jacobian for the gas-problem...
std::vector< size_t > m_spSurfLarge
Vector containing the indices of the largest species in each surface phase.
vector_fp m_CSolnSave
Temporary vector with length equal to max m_maxTotSpecies.
std::vector< size_t > m_nSpeciesSurfPhase
Vector of length number of surface phases containing the number of surface species in each phase...
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
void updateMFKinSpecies(doublereal *XMolKinSp, int isp)
Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceK...
Base class for exceptions thrown by Cantera classes.
vector_fp m_numEqn2
Temporary vector with length equal to max m_maxTotSpecies.
const int SFLUX_TRANSIENT
The transient calculation is performed here for an amount of time specified by "time_scale".
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
std::vector< SurfPhase * > m_ptrsSurfPhase
Vector of surface phase pointers.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
std::vector< size_t > m_kinObjPhaseIDSurfPhase
Phase ID in the InterfaceKinetics object of the surface phase.
vector_fp m_netProductionRatesSave
Temporary vector with length equal to max m_maxTotSpecies.
std::vector< size_t > m_kinObjIndex
Index between the equation index and the index of the InterfaceKinetics object.
Advances the surface coverages of the associated set of SurfacePhase objects in time.
size_t m_numSurfPhases
Number of surface phases in the surface problem.
std::vector< ThermoPhase * > m_bulkPhasePtrs
Vector of bulk phase pointers, length is equal to m_numBulkPhases.
std::vector< size_t > m_numBulkSpecies
Vector of number of species in the m_numBulkPhases phases.
vector_fp m_CSolnSPOld
Saved solution vector at the old time step. length MAX(1, m_neq)
vector_fp m_resid
Residual for the surface problem.
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).
std::vector< size_t > m_kinSpecIndex
Index between the equation index and the position in the kinetic species array for the appropriate ki...
size_t m_maxTotSpecies
Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution s...
DenseMatrix m_Jac
Jacobian.
std::vector< size_t > m_indexKinObjSurfPhase
Mapping between the surface phases and the InterfaceKinetics objects.
vector_fp m_wtSpecies
Weights for the species concentrations norm calculation.
Namespace for the Cantera kernel.
void fun_eval(doublereal *resid, const doublereal *CSolnSP, const doublereal *CSolnOldSP, const bool do_time, const doublereal deltaT)
Main Function evaluation.
size_t m_numTotBulkSpeciesSS
Total number of species in all bulk phases.
void updateState(const doublereal *cSurfSpec)
Update the surface states of the surface phases.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
void calcWeights(doublereal wtSpecies[], doublereal wtResid[], const Array2D &Jac, const doublereal CSolnSP[], const doublereal abstol, const doublereal reltol)
Calculate the solution and residual weights.