Cantera  3.1.0a1
Loading...
Searching...
No Matches
solveSP.cpp
Go to the documentation of this file.
1/**
2 * @file: solveSP.cpp Implicit surface site concentration solver
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
11
12namespace Cantera
13{
14
15// STATIC ROUTINES DEFINED IN THIS FILE
16
17static double calc_damping(double* x, double* dx, size_t dim, int*);
18static double calcWeightedNorm(const double [], const double dx[], size_t);
19
20// solveSP Class Definitions
21
22solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
23 m_objects(surfChemPtr->getObjects()),
24 m_bulkFunc(bulkFunc)
25{
26 for (size_t n = 0; n < m_objects.size(); n++) {
28 SurfPhase* sp = dynamic_cast<SurfPhase*>(&kin->thermo(0));
29 if (sp == nullptr) {
30 throw CanteraError("solveSP::solveSP",
31 "InterfaceKinetics object has no surface phase");
32 }
33
35 m_indexKinObjSurfPhase.push_back(n);
36
37 m_ptrsSurfPhase.push_back(sp);
38 size_t nsp = sp->nSpecies();
39 m_nSpeciesSurfPhase.push_back(nsp);
41 }
42
43 if (bulkFunc == BULK_DEPOSITION) {
45 } else {
47 }
48
49 for (size_t n = 0; n < m_numSurfPhases; n++) {
50 size_t tsp = m_objects[n]->nTotalSpecies();
51 m_maxTotSpecies = std::max(m_maxTotSpecies, tsp);
52 }
54
56 m_numEqn1.resize(m_maxTotSpecies, 0.0);
57 m_numEqn2.resize(m_maxTotSpecies, 0.0);
59 m_CSolnSave.resize(m_neq, 0.0);
64
65 size_t kindexSP = 0;
66 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
67 size_t iKinObject = m_indexKinObjSurfPhase[isp];
68 InterfaceKinetics* kin = m_objects[iKinObject];
69 size_t kstart = kin->kineticsSpeciesIndex(0, 0);
70 size_t nsp = m_nSpeciesSurfPhase[isp];
71 m_eqnIndexStartSolnPhase[isp] = kindexSP;
72 for (size_t k = 0; k < nsp; k++, kindexSP++) {
73 m_kinSpecIndex[kindexSP] = kstart + k;
74 m_kinObjIndex[kindexSP] = isp;
75 }
76 }
77
78 // Dimension solution vector
79 size_t dim1 = std::max<size_t>(1, m_neq);
80 m_CSolnSP.resize(dim1, 0.0);
81 m_CSolnSPInit.resize(dim1, 0.0);
82 m_CSolnSPOld.resize(dim1, 0.0);
83 m_wtResid.resize(dim1, 0.0);
84 m_wtSpecies.resize(dim1, 0.0);
85 m_resid.resize(dim1, 0.0);
86 m_Jac.resize(dim1, dim1, 0.0);
87}
88
89int solveSP::solveSurfProb(int ifunc, double time_scale, double TKelvin,
90 double PGas, double reltol, double abstol)
91{
92 double EXTRA_ACCURACY = 0.001;
93 if (ifunc == SFLUX_JACOBIAN) {
94 EXTRA_ACCURACY *= 0.001;
95 }
96 int label_t=-1; // Species IDs for time control
97 int label_d = -1; // Species IDs for damping control
98 int label_t_old=-1;
99 double label_factor = 1.0;
100 int iter=0; // iteration number on nonlinear solver
101 int iter_max=1000; // maximum number of nonlinear iterations
102 double deltaT = 1.0E-10; // Delta time step
103 double damp=1.0;
104 double inv_t = 0.0;
105 double t_real = 0.0, update_norm = 1.0E6;
106 bool do_time = false, not_converged = true;
107 m_ioflag = std::min(m_ioflag, 1);
108
109 // Set the initial value of the do_time parameter
110 if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) {
111 do_time = true;
112 }
113
114 // Store the initial guess for the surface problem in the soln vector,
115 // CSoln, and in an separate vector CSolnInit.
116 size_t loc = 0;
117 for (size_t n = 0; n < m_numSurfPhases; n++) {
118 m_ptrsSurfPhase[n]->getConcentrations(m_numEqn1.data());
119 for (size_t k = 0; k < m_nSpeciesSurfPhase[n]; k++) {
120 m_CSolnSP[loc] = m_numEqn1[k];
121 loc++;
122 }
123 }
124
126
127 // Calculate the largest species in each phase
128 evalSurfLarge(m_CSolnSP.data());
129
130 if (m_ioflag) {
131 print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol);
132 }
133
134 // Quick return when there isn't a surface problem to solve
135 if (m_neq == 0) {
136 not_converged = false;
137 update_norm = 0.0;
138 }
139
140 // Start of Newton's method
141 while (not_converged && iter < iter_max) {
142 iter++;
143 // Store previous iteration's solution in the old solution vector
145
146 // Evaluate the largest surface species for each surface phase every
147 // 5 iterations.
148 if (iter%5 == 4) {
149 evalSurfLarge(m_CSolnSP.data());
150 }
151
152 // Calculate the value of the time step
153 // - heuristics to stop large oscillations in deltaT
154 if (do_time) {
155 // don't hurry increase in time step at the same time as damping
156 if (damp < 1.0) {
157 label_factor = 1.0;
158 }
159 double tmp = calc_t(m_netProductionRatesSave.data(),
160 m_XMolKinSpecies.data(),
161 &label_t, &label_t_old, &label_factor, m_ioflag);
162 if (iter < 10) {
163 inv_t = tmp;
164 } else if (tmp > 2.0*inv_t) {
165 inv_t = 2.0*inv_t;
166 } else {
167 inv_t = tmp;
168 }
169
170 // Check end condition
171 if (ifunc == SFLUX_TRANSIENT) {
172 tmp = t_real + 1.0/inv_t;
173 if (tmp > time_scale) {
174 inv_t = 1.0/(time_scale - t_real);
175 }
176 }
177 } else {
178 // make steady state calc a step of 1 million seconds to prevent
179 // singular Jacobians for some pathological cases
180 inv_t = 1.0e-6;
181 }
182 deltaT = 1.0/inv_t;
183
184 // Call the routine to numerically evaluation the Jacobian and residual
185 // for the current iteration.
186 resjac_eval(m_Jac, m_resid.data(), m_CSolnSP.data(),
187 m_CSolnSPOld.data(), do_time, deltaT);
188
189 // Calculate the weights. Make sure the calculation is carried out on
190 // the first iteration.
191 if (iter%4 == 1) {
192 calcWeights(m_wtSpecies.data(), m_wtResid.data(),
193 m_Jac, m_CSolnSP.data(), abstol, reltol);
194 }
195
196 // Find the weighted norm of the residual
197 double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
198
199 // Solve Linear system. The solution is in m_resid
200 solve(m_Jac, m_resid.data());
201
202 // Calculate the Damping factor needed to keep all unknowns between 0
203 // and 1, and not allow too large a change (factor of 2) in any unknown.
204 damp = calc_damping(m_CSolnSP.data(), m_resid.data(), m_neq, &label_d);
205
206 // Calculate the weighted norm of the update vector Here, resid is the
207 // delta of the solution, in concentration units.
208 update_norm = calcWeightedNorm(m_wtSpecies.data(),
209 m_resid.data(), m_neq);
210
211 // Update the solution vector and real time Crop the concentrations to
212 // zero.
213 for (size_t irow = 0; irow < m_neq; irow++) {
214 m_CSolnSP[irow] -= damp * m_resid[irow];
215 }
216 for (size_t irow = 0; irow < m_neq; irow++) {
217 m_CSolnSP[irow] = std::max(0.0, m_CSolnSP[irow]);
218 }
219 updateState(m_CSolnSP.data());
220
221 if (do_time) {
222 t_real += damp/inv_t;
223 }
224
225 if (m_ioflag) {
226 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
227 update_norm, resid_norm, do_time);
228 }
229
230 if (ifunc == SFLUX_TRANSIENT) {
231 not_converged = (t_real < time_scale);
232 } else {
233 if (do_time) {
234 if (t_real > time_scale ||
235 (resid_norm < 1.0e-7 &&
236 update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
237 do_time = false;
238 }
239 } else {
240 not_converged = ((update_norm > EXTRA_ACCURACY) ||
241 (resid_norm > EXTRA_ACCURACY));
242 }
243 }
244 } // End of Newton's Method while statement
245
246 // End Newton's method. If not converged, print error message and
247 // recalculate sdot's at equal site fractions.
248 if (not_converged && m_ioflag) {
249 writelog("#$#$#$# Error in solveSP $#$#$#$ \n");
250 writelogf("Newton iter on surface species did not converge, "
251 "update_norm = %e \n", update_norm);
252 writelog("Continuing anyway\n");
253 }
254
255 // Decide on what to return in the solution vector. Right now, will always
256 // return the last solution no matter how bad
257 if (m_ioflag) {
258 fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(),
259 false, deltaT);
260 double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
261 printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
262 update_norm, resid_norm, do_time, true);
263 }
264
265 // Return with the appropriate flag
266 if (update_norm > 1.0) {
267 return -1;
268 }
269 return 1;
270}
271
272void solveSP::updateState(const double* CSolnSP)
273{
274 vector<double> X;
275 size_t loc = 0;
276 for (size_t n = 0; n < m_numSurfPhases; n++) {
277 X.resize(m_nSpeciesSurfPhase[n]);
278 for (size_t k = 0; k < X.size(); k++) {
279 X[k] = CSolnSP[loc + k] / m_ptrsSurfPhase[n]->siteDensity();
280 }
281 m_ptrsSurfPhase[n]->setMoleFractions_NoNorm(X.data());
282 loc += m_nSpeciesSurfPhase[n];
283 }
284}
285
286void solveSP::updateMFSolnSP(double* XMolSolnSP)
287{
288 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
289 size_t keqnStart = m_eqnIndexStartSolnPhase[isp];
290 m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart);
291 }
292}
293
294void solveSP::updateMFKinSpecies(double* XMolKinSpecies, int isp)
295{
296 InterfaceKinetics* kin = m_objects[isp];
297 for (size_t iph = 0; iph < kin->nPhases(); iph++) {
298 size_t ksi = kin->kineticsSpeciesIndex(0, iph);
299 kin->thermo(iph).getMoleFractions(XMolKinSpecies + ksi);
300 }
301}
302
303void solveSP::evalSurfLarge(const double* CSolnSP)
304{
305 size_t kindexSP = 0;
306 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
307 double Clarge = CSolnSP[kindexSP];
308 m_spSurfLarge[isp] = 0;
309 kindexSP++;
310 for (size_t k = 1; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
311 if (CSolnSP[kindexSP] > Clarge) {
312 Clarge = CSolnSP[kindexSP];
313 m_spSurfLarge[isp] = k;
314 }
315 }
316 }
317}
318
319void solveSP::fun_eval(double* resid, const double* CSoln, const double* CSolnOld,
320 const bool do_time, const double deltaT)
321{
322 size_t k;
323 double lenScale = 1.0E-9;
324 if (m_numSurfPhases > 0) {
325 // update the surface concentrations with the input surface
326 // concentration vector
327 updateState(CSoln);
328
329 // Get the net production rates of all of the species in the
330 // surface kinetics mechanism
331 //
332 // HKM Should do it here for all kinetics objects so that
333 // bulk will eventually work.
334 if (do_time) {
335 size_t kindexSP = 0;
336 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
337 size_t nsp = m_nSpeciesSurfPhase[isp];
338 InterfaceKinetics* kinPtr = m_objects[isp];
339 size_t kins = kindexSP;
341 for (k = 0; k < nsp; k++, kindexSP++) {
342 resid[kindexSP] =
343 (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
345 }
346
347 size_t kspecial = kins + m_spSurfLarge[isp];
348 double sd = m_ptrsSurfPhase[isp]->siteDensity();
349 resid[kspecial] = sd;
350 for (k = 0; k < nsp; k++) {
351 resid[kspecial] -= CSoln[kins + k];
352 }
353 }
354 } else {
355 size_t kindexSP = 0;
356 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
357 size_t nsp = m_nSpeciesSurfPhase[isp];
358 InterfaceKinetics* kinPtr = m_objects[isp];
359 size_t kins = kindexSP;
361 for (k = 0; k < nsp; k++, kindexSP++) {
362 resid[kindexSP] = - m_netProductionRatesSave[k];
363 }
364 size_t kspecial = kins + m_spSurfLarge[isp];
365 double sd = m_ptrsSurfPhase[isp]->siteDensity();
366 resid[kspecial] = sd;
367 for (k = 0; k < nsp; k++) {
368 resid[kspecial] -= CSoln[kins + k];
369 }
370 }
371 }
372
374 size_t kindexSP = m_numTotSurfSpecies;
375 for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
376 double* XBlk = m_numEqn1.data();
377 size_t nsp = m_nSpeciesSurfPhase[isp];
378 double grRate = 0.0;
379 for (k = 0; k < nsp; k++) {
380 if (m_netProductionRatesSave[k] > 0.0) {
381 grRate += m_netProductionRatesSave[k];
382 }
383 }
384 resid[kindexSP] = m_bulkPhasePtrs[isp]->molarDensity();
385 for (k = 0; k < nsp; k++) {
386 resid[kindexSP] -= CSoln[kindexSP + k];
387 }
388 if (grRate > 0.0) {
389 for (k = 1; k < nsp; k++) {
390 if (m_netProductionRatesSave[k] > 0.0) {
391 resid[kindexSP + k] = XBlk[k] * grRate
393 } else {
394 resid[kindexSP + k] = XBlk[k] * grRate;
395 }
396 }
397 } else {
398 grRate = 1.0E-6;
399 //! @todo the appearance of k in this formula is suspicious
400 grRate += fabs(m_netProductionRatesSave[k]);
401 for (k = 1; k < nsp; k++) {
402 resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp);
403 }
404 }
405 if (do_time) {
406 for (k = 1; k < nsp; k++) {
407 resid[kindexSP + k] +=
408 lenScale / deltaT *
409 (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]);
410 }
411 }
412 kindexSP += nsp;
413 }
414 }
415 }
416}
417
418void solveSP::resjac_eval(DenseMatrix& jac, double resid[], double CSoln[],
419 const double CSolnOld[], const bool do_time,
420 const double deltaT)
421{
422 size_t kColIndex = 0;
423 // Calculate the residual
424 fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
425 // Now we will look over the columns perturbing each unknown.
426 for (size_t jsp = 0; jsp < m_numSurfPhases; jsp++) {
427 size_t nsp = m_nSpeciesSurfPhase[jsp];
428 double sd = m_ptrsSurfPhase[jsp]->siteDensity();
429 for (size_t kCol = 0; kCol < nsp; kCol++) {
430 double cSave = CSoln[kColIndex];
431 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
432 CSoln[kColIndex] += dc;
433 fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
434 for (size_t i = 0; i < m_neq; i++) {
435 jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
436 }
437 CSoln[kColIndex] = cSave;
438 kColIndex++;
439 }
440 }
441
443 for (size_t jsp = 0; jsp < m_numBulkPhasesSS; jsp++) {
444 size_t nsp = m_numBulkSpecies[jsp];
445 double sd = m_bulkPhasePtrs[jsp]->molarDensity();
446 for (size_t kCol = 0; kCol < nsp; kCol++) {
447 double cSave = CSoln[kColIndex];
448 double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
449 CSoln[kColIndex] += dc;
450 fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
451 for (size_t i = 0; i < m_neq; i++) {
452 jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
453 }
454 CSoln[kColIndex] = cSave;
455 kColIndex++;
456 }
457 }
458 }
459}
460
461/**
462 * This function calculates a damping factor for the Newton iteration update
463 * vector, dxneg, to insure that all site and bulk fractions, x, remain
464 * bounded between zero and one.
465 *
466 * dxneg[] = negative of the update vector.
467 *
468 * The constant "APPROACH" sets the fraction of the distance to the boundary
469 * that the step can take. If the full step would not force any fraction
470 * outside of 0-1, then Newton's method is allowed to operate normally.
471 */
472static double calc_damping(double x[], double dxneg[], size_t dim, int* label)
473{
474 const double APPROACH = 0.80;
475 double damp = 1.0;
476 static double damp_old = 1.0; //! @todo this variable breaks thread safety
477 *label = -1;
478
479 for (size_t i = 0; i < dim; i++) {
480 // Calculate the new suggested new value of x[i]
481 double xnew = x[i] - damp * dxneg[i];
482
483 // Calculate the allowed maximum and minimum values of x[i]
484 // - Only going to allow x[i] to converge to zero by a
485 // single order of magnitude at a time
486 double xtop = 1.0 - 0.1*fabs(1.0-x[i]);
487 double xbot = fabs(x[i]*0.1) - 1.0e-16;
488 if (xnew > xtop) {
489 damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
490 *label = int(i);
491 } else if (xnew < xbot) {
492 damp = APPROACH * x[i] / dxneg[i];
493 *label = int(i);
494 } else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
495 damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
496 *label = int(i);
497 }
498 }
499 damp = std::max(damp, 1e-2);
500
501 // Only allow the damping parameter to increase by a factor of three each
502 // iteration. Heuristic to avoid oscillations in the value of damp
503 if (damp > damp_old*3) {
504 damp = damp_old*3;
505 *label = -1;
506 }
507
508 // Save old value of the damping parameter for use in subsequent calls.
509 damp_old = damp;
510 return damp;
511
512} /* calc_damping */
513
514/**
515 * This function calculates the norm of an update, dx[], based on the
516 * weighted values of x.
517 */
518static double calcWeightedNorm(const double wtX[], const double dx[], size_t dim)
519{
520 double norm = 0.0;
521 if (dim == 0) {
522 return 0.0;
523 }
524 for (size_t i = 0; i < dim; i++) {
525 norm += pow(dx[i] / wtX[i], 2);
526 }
527 return sqrt(norm/dim);
528}
529
530void solveSP::calcWeights(double wtSpecies[], double wtResid[],
531 const Array2D& Jac, const double CSoln[],
532 const double abstol, const double reltol)
533{
534 // First calculate the weighting factor for the concentrations of the
535 // surface species and bulk species.
536 size_t kindex = 0;
537 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
538 double sd = m_ptrsSurfPhase[isp]->siteDensity();
539 for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindex++) {
540 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
541 }
542 }
544 for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
545 double sd = m_bulkPhasePtrs[isp]->molarDensity();
546 for (size_t k = 0; k < m_numBulkSpecies[isp]; k++, kindex++) {
547 wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
548 }
549 }
550 }
551
552 // Now do the residual Weights. Since we have the Jacobian, we will use it
553 // to generate a number based on the what a significant change in a solution
554 // variable does to each residual. This is a row sum scale operation.
555 for (size_t k = 0; k < m_neq; k++) {
556 wtResid[k] = 0.0;
557 for (size_t jcol = 0; jcol < m_neq; jcol++) {
558 wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
559 }
560 }
561}
562
563double solveSP::calc_t(double netProdRateSolnSP[], double XMolSolnSP[], int* label,
564 int* label_old, double* label_factor, int ioflag)
565{
566 double inv_timeScale = 1.0E-10;
567 size_t kindexSP = 0;
568 *label = 0;
569 updateMFSolnSP(XMolSolnSP);
570 for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
571 // Get the interface kinetics associated with this surface
572 InterfaceKinetics* kin = m_objects[isp];
573
574 kin->getNetProductionRates(m_numEqn1.data());
575 double sden = kin->thermo(0).molarDensity();
576 for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
577 netProdRateSolnSP[kindexSP] = m_numEqn1[k];
578 double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10);
579 tmp *= sden;
580 tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
581 if (netProdRateSolnSP[kindexSP]> 0.0) {
582 tmp /= 100.;
583 }
584 if (tmp > inv_timeScale) {
585 inv_timeScale = tmp;
586 *label = int(kindexSP);
587 }
588 }
589 }
590
591 // Increase time step exponentially as same species repeatedly controls time
592 // step
593 if (*label == *label_old) {
594 *label_factor *= 1.5;
595 } else {
596 *label_old = *label;
597 *label_factor = 1.0;
598 }
599 return inv_timeScale / *label_factor;
600} // calc_t
601
602void solveSP::print_header(int ioflag, int ifunc, double time_scale,
603 int damping, double reltol, double abstol)
604{
605 if (ioflag) {
606 writelog("\n================================ SOLVESP CALL SETUP "
607 "========================================\n");
608 if (ifunc == SFLUX_INITIALIZE) {
609 writelog("\n SOLVESP Called with Initialization turned on\n");
610 writelogf(" Time scale input = %9.3e\n", time_scale);
611 } else if (ifunc == SFLUX_RESIDUAL) {
612 writelog("\n SOLVESP Called to calculate steady state residual\n");
613 writelog(" from a good initial guess\n");
614 } else if (ifunc == SFLUX_JACOBIAN) {
615 writelog("\n SOLVESP Called to calculate steady state Jacobian\n");
616 writelog(" from a good initial guess\n");
617 } else if (ifunc == SFLUX_TRANSIENT) {
618 writelog("\n SOLVESP Called to integrate surface in time\n");
619 writelogf(" for a total of %9.3e sec\n", time_scale);
620 } else {
621 throw CanteraError("solveSP::print_header",
622 "Unknown ifunc flag = {}", ifunc);
623 }
624
626 writelog(" The composition of the Bulk Phases will be calculated\n");
627 } else if (m_bulkFunc == BULK_ETCH) {
628 writelog(" Bulk Phases have fixed compositions\n");
629 } else {
630 throw CanteraError("solveSP::print_header",
631 "Unknown bulkFunc flag = {}", m_bulkFunc);
632 }
633
634 if (damping) {
635 writelog(" Damping is ON \n");
636 } else {
637 writelog(" Damping is OFF \n");
638 }
639
640 writelogf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
641 }
642
643 if (ioflag == 1) {
644 writelog("\n\n\t Iter Time Del_t Damp DelX "
645 " Resid Name-Time Name-Damp\n");
646 writelog("\t -----------------------------------------------"
647 "------------------------------------\n");
648 }
649}
650
651void solveSP::printIteration(int ioflag, double damp, int label_d,
652 int label_t, double inv_t, double t_real,
653 size_t iter, double update_norm,
654 double resid_norm, bool do_time, bool final)
655{
656 if (ioflag == 1) {
657 if (final) {
658 writelogf("\tFIN%3d ", iter);
659 } else {
660 writelogf("\t%6d ", iter);
661 }
662 if (do_time) {
663 writelogf("%9.4e %9.4e ", t_real, 1.0/inv_t);
664 } else {
665 writeline(' ', 22, false);
666 }
667 if (damp < 1.0) {
668 writelogf("%9.4e ", damp);
669 } else {
670 writeline(' ', 11, false);
671 }
672 writelogf("%9.4e %9.4e", update_norm, resid_norm);
673 if (do_time) {
674 size_t k = m_kinSpecIndex[label_t];
675 size_t isp = m_kinObjIndex[label_t];
676 writelog(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
677 } else {
678 writeline(' ', 16, false);
679 }
680 if (label_d >= 0) {
681 size_t k = m_kinSpecIndex[label_d];
682 size_t isp = m_kinObjIndex[label_d];
683 writelogf(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
684 }
685 if (final) {
686 writelog(" -- success");
687 }
688 writelog("\n");
689 }
690} // printIteration
691
692}
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.
Definition Array.h:32
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
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.
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition Kinetics.h:242
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition Kinetics.h:184
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition Kinetics.h:276
virtual void getNetProductionRates(double *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition Kinetics.cpp:363
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:576
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:231
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:434
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition SurfPhase.h:98
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
solveSP(ImplicitSurfChem *surfChemPtr, int bulkFunc=BULK_ETCH)
Constructor for the object.
Definition solveSP.cpp:22
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
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
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
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
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.
Definition solveSP.cpp:518
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).