17static double calc_damping(
double* x,
double* dx,
size_t dim,
int*);
23 m_objects(surfChemPtr->getObjects()),
26 for (
size_t n = 0; n <
m_objects.size(); n++) {
32 "InterfaceKinetics object has no surface phase");
52 size_t tsp =
m_objects[n]->nTotalSpecies();
75 for (
size_t k = 0; k < nsp; k++, kindexSP++) {
82 size_t dim1 = std::max<size_t>(1,
m_neq);
93 double PGas,
double reltol,
double abstol)
95 double EXTRA_ACCURACY = 0.001;
97 EXTRA_ACCURACY *= 0.001;
102 double label_factor = 1.0;
105 double deltaT = 1.0E-10;
108 double t_real = 0.0, update_norm = 1.0E6;
109 bool do_time =
false, not_converged =
true;
110 m_ioflag = std::min(m_ioflag, 1);
134 print_header(m_ioflag, ifunc, time_scale,
true, reltol, abstol);
139 not_converged =
false;
144 while (not_converged && iter < iter_max) {
164 &label_t, &label_t_old, &label_factor, m_ioflag);
167 }
else if (tmp > 2.0*inv_t) {
175 tmp = t_real + 1.0/inv_t;
176 if (tmp > time_scale) {
177 inv_t = 1.0/(time_scale - t_real);
216 for (
size_t irow = 0; irow <
m_neq; irow++) {
219 for (
size_t irow = 0; irow <
m_neq; irow++) {
225 t_real += damp/inv_t;
229 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
230 update_norm, resid_norm, do_time);
234 not_converged = (t_real < time_scale);
237 if (t_real > time_scale ||
238 (resid_norm < 1.0e-7 &&
239 update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
243 not_converged = ((update_norm > EXTRA_ACCURACY) ||
244 (resid_norm > EXTRA_ACCURACY));
251 if (not_converged && m_ioflag) {
252 writelog(
"#$#$#$# Error in solveSP $#$#$#$ \n");
253 writelogf(
"Newton iter on surface species did not converge, "
254 "update_norm = %e \n", update_norm);
264 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
265 update_norm, resid_norm, do_time,
true);
269 if (update_norm > 1.0) {
281 for (
size_t k = 0; k < X.size(); k++) {
300 for (
size_t iph = 0; iph < kin->
nPhases(); iph++) {
310 double Clarge = CSolnSP[kindexSP];
314 if (CSolnSP[kindexSP] > Clarge) {
315 Clarge = CSolnSP[kindexSP];
323 const bool do_time,
const double deltaT)
326 double lenScale = 1.0E-9;
344 size_t kins = kindexSP;
346 for (k = 0; k < nsp; k++, kindexSP++) {
348 (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
354 resid[kspecial] = sd;
355 for (k = 0; k < nsp; k++) {
356 resid[kspecial] -= CSoln[kins + k];
366 size_t kins = kindexSP;
368 for (k = 0; k < nsp; k++, kindexSP++) {
373 resid[kspecial] = sd;
374 for (k = 0; k < nsp; k++) {
375 resid[kspecial] -= CSoln[kins + k];
389 for (k = 0; k < nsp; k++) {
395 for (k = 0; k < nsp; k++) {
396 resid[kindexSP] -= CSoln[kindexSP + k];
399 for (k = 1; k < nsp; k++) {
401 resid[kindexSP + k] = XBlk[k] * grRate
404 resid[kindexSP + k] = XBlk[k] * grRate;
411 for (k = 1; k < nsp; k++) {
412 resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp);
416 for (k = 1; k < nsp; k++) {
417 resid[kindexSP + k] +=
419 (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]);
429 const double CSolnOld[],
const bool do_time,
432 size_t kColIndex = 0;
434 fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
439 for (
size_t kCol = 0; kCol < nsp; kCol++) {
440 double cSave = CSoln[kColIndex];
441 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
442 CSoln[kColIndex] += dc;
444 for (
size_t i = 0; i <
m_neq; i++) {
445 jac(i, kColIndex) = (
m_numEqn2[i] - resid[i])/dc;
447 CSoln[kColIndex] = cSave;
456 for (
size_t kCol = 0; kCol < nsp; kCol++) {
457 double cSave = CSoln[kColIndex];
458 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
459 CSoln[kColIndex] += dc;
461 for (
size_t i = 0; i <
m_neq; i++) {
462 jac(i, kColIndex) = (
m_numEqn2[i] - resid[i])/dc;
464 CSoln[kColIndex] = cSave;
482static double calc_damping(
double x[],
double dxneg[],
size_t dim,
int* label)
484 const double APPROACH = 0.80;
486 static double damp_old = 1.0;
489 for (
size_t i = 0; i < dim; i++) {
491 double xnew = x[i] - damp * dxneg[i];
496 double xtop = 1.0 - 0.1*fabs(1.0-x[i]);
497 double xbot = fabs(x[i]*0.1) - 1.0e-16;
499 damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
501 }
else if (xnew < xbot) {
502 damp = APPROACH * x[i] / dxneg[i];
504 }
else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
505 damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
509 damp = std::max(damp, 1e-2);
513 if (damp > damp_old*3) {
534 for (
size_t i = 0; i < dim; i++) {
535 norm += pow(dx[i] / wtX[i], 2);
537 return sqrt(norm/dim);
541 const Array2D& Jac,
const double CSoln[],
542 const double abstol,
const double reltol)
550 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
557 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
565 for (
size_t k = 0; k <
m_neq; k++) {
567 for (
size_t jcol = 0; jcol <
m_neq; jcol++) {
568 wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
574 int* label_old,
double* label_factor,
int ioflag)
576 double inv_timeScale = 1.0E-10;
591 size_t kspindex = kstart + k;
592 netProdRateSolnSP[kindexSP] =
m_numEqn1[kspindex];
593 double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10);
595 tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
596 if (netProdRateSolnSP[kindexSP]> 0.0) {
599 if (tmp > inv_timeScale) {
601 *label = int(kindexSP);
608 if (*label == *label_old) {
609 *label_factor *= 1.5;
614 return inv_timeScale / *label_factor;
618 int damping,
double reltol,
double abstol)
621 writelog(
"\n================================ SOLVESP CALL SETUP "
622 "========================================\n");
624 writelog(
"\n SOLVESP Called with Initialization turned on\n");
625 writelogf(
" Time scale input = %9.3e\n", time_scale);
627 writelog(
"\n SOLVESP Called to calculate steady state residual\n");
628 writelog(
" from a good initial guess\n");
630 writelog(
"\n SOLVESP Called to calculate steady state Jacobian\n");
631 writelog(
" from a good initial guess\n");
633 writelog(
"\n SOLVESP Called to integrate surface in time\n");
634 writelogf(
" for a total of %9.3e sec\n", time_scale);
637 "Unknown ifunc flag = {}", ifunc);
641 writelog(
" The composition of the Bulk Phases will be calculated\n");
643 writelog(
" Bulk Phases have fixed compositions\n");
655 writelogf(
" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
659 writelog(
"\n\n\t Iter Time Del_t Damp DelX "
660 " Resid Name-Time Name-Damp\n");
661 writelog(
"\t -----------------------------------------------"
662 "------------------------------------\n");
667 int label_t,
double inv_t,
double t_real,
668 size_t iter,
double update_norm,
669 double resid_norm,
bool do_time,
bool final)
678 writelogf(
"%9.4e %9.4e ", t_real, 1.0/inv_t);
680 writeline(
' ', 22,
false);
685 writeline(
' ', 11,
false);
687 writelogf(
"%9.4e %9.4e", update_norm, resid_norm);
693 writeline(
' ', 16,
false);
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
Advances the surface coverages of the associated set of SurfacePhase objects in time.
A kinetics manager for heterogeneous reaction mechanisms.
size_t reactionPhaseIndex() const
Phase where the reactions occur.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
virtual double molarDensity() const
Molar density (kmol/m^3).
size_t nSpecies() const
Returns the number of species in the phase.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
vector< size_t > m_numBulkSpecies
Vector of number of species in the m_numBulkPhases phases.
vector< double > m_wtSpecies
Weights for the species concentrations norm calculation.
size_t m_numTotSurfSpecies
Total number of surface species in all surface phases.
void evalSurfLarge(const double *CSolnSP)
Update the vector that keeps track of the largest species in each surface phase.
void updateState(const double *cSurfSpec)
Update the surface states of the surface phases.
vector< size_t > m_indexKinObjSurfPhase
Mapping between the surface phases and the InterfaceKinetics objects.
vector< size_t > m_kinObjIndex
Index between the equation index and the index of the InterfaceKinetics object.
solveSP(ImplicitSurfChem *surfChemPtr, int bulkFunc=BULK_ETCH)
Constructor for the object.
vector< double > m_resid
Residual for the surface problem.
vector< ThermoPhase * > m_bulkPhasePtrs
Vector of bulk phase pointers, length is equal to m_numBulkPhases.
vector< size_t > m_spSurfLarge
Vector containing the indices of the largest species in each surface phase.
vector< double > m_netProductionRatesSave
Temporary vector with length equal to max m_maxTotSpecies.
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.
int m_bulkFunc
This variable determines how the bulk phases are to be handled.
vector< double > m_wtResid
Weights for the residual norm calculation. length MAX(1, m_neq)
vector< InterfaceKinetics * > & m_objects
Vector of interface kinetics objects.
void fun_eval(double *resid, const double *CSolnSP, const double *CSolnOldSP, const bool do_time, const double deltaT)
Main Function evaluation.
vector< double > m_CSolnSave
Temporary vector with length equal to max m_maxTotSpecies.
void updateMFSolnSP(double *XMolSolnSP)
Update mole fraction vector consisting of unknowns in surface problem.
vector< double > m_CSolnSPOld
Saved solution vector at the old time step. length MAX(1, m_neq)
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.
vector< double > m_CSolnSP
Solution vector. length MAX(1, m_neq)
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_numTotBulkSpeciesSS
Total number of species in all bulk phases.
vector< double > m_numEqn1
Temporary vector with length equal to max m_maxTotSpecies.
vector< double > m_CSolnSPInit
Saved initial solution vector. length MAX(1, m_neq)
vector< double > m_XMolKinSpecies
Vector of mole fractions. length m_maxTotSpecies.
vector< double > m_numEqn2
Temporary vector with length equal to max m_maxTotSpecies.
vector< size_t > m_eqnIndexStartSolnPhase
Index of the start of the unknowns for each solution phase.
vector< size_t > m_nSpeciesSurfPhase
Vector of length number of surface phases containing the number of surface species in each phase.
size_t m_numSurfPhases
Number of surface phases in the surface problem.
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.
void updateMFKinSpecies(double *XMolKinSp, int isp)
Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceK...
DenseMatrix m_Jac
Jacobian.
size_t m_numBulkPhasesSS
Total number of volumetric condensed phases included in the steady state problem handled by this rout...
size_t m_neq
Total number of equations to solve in the implicit problem.
vector< SurfPhase * > m_ptrsSurfPhase
Vector of surface phase pointers.
vector< size_t > m_kinObjPhaseIDSurfPhase
Phase ID in the InterfaceKinetics object of the surface phase.
size_t m_maxTotSpecies
Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution s...
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.
void calcWeights(double wtSpecies[], double wtResid[], const Array2D &Jac, const double CSolnSP[], const double abstol, const double reltol)
Calculate the solution and residual weights.
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.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
const int BULK_ETCH
Etching of a bulk phase is to be expected.
const int BULK_DEPOSITION
Deposition of a bulk phase is to be expected.
const int SFLUX_TRANSIENT
The transient calculation is performed here for an amount of time specified by "time_scale".
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species.
const int SFLUX_JACOBIAN
Calculation of the surface problem is due to the need for a numerical Jacobian for the gas-problem.
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one.
Namespace for the Cantera kernel.
static double calcWeightedNorm(const double[], const double dx[], size_t)
This function calculates the norm of an update, dx[], based on the weighted values of x.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).