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