Cantera  2.0
InterfaceKinetics.cpp
Go to the documentation of this file.
1 /**
2  * @file InterfaceKinetics.cpp
3  *
4  */
5 
6 // Copyright 2002 California Institute of Technology
7 
11 
14 
15 #include "ImplicitSurfChem.h"
16 
17 using namespace std;
18 
19 namespace Cantera
20 {
21 //====================================================================================================================
22 InterfaceKineticsData::InterfaceKineticsData() :
23  m_logp0(0.0),
24  m_logc0(0.0),
25  m_ROP_ok(false),
26  m_temp(0.0),
27  m_logtemp(0.0)
28 {
29 }
30 //====================================================================================================================
31 InterfaceKineticsData:: InterfaceKineticsData(const InterfaceKineticsData& right) :
32  m_logp0(0.0),
33  m_logc0(0.0),
34  m_ROP_ok(false),
35  m_temp(0.0),
36  m_logtemp(0.0)
37 {
38  *this = right;
39 }
40 //====================================================================================================================
41 InterfaceKineticsData::~InterfaceKineticsData()
42 {
43 }
44 //====================================================================================================================
45 InterfaceKineticsData& InterfaceKineticsData::operator=(const InterfaceKineticsData& right)
46 {
47  if (this == &right) {
48  return *this;
49  }
50  m_logp0 = right.m_logp0;
51  m_logc0 = right.m_logc0;
52  m_ropf = right.m_ropf;
53  m_ropr = right.m_ropr;
54  m_ropnet = right.m_ropnet;
55  m_ROP_ok = right.m_ROP_ok;
56  m_temp = right.m_temp;
57  m_logtemp = right.m_logtemp;
58  m_rfn = right.m_rfn;
59  m_rkcn = right.m_rkcn;
60  return *this;
61 }
62 //====================================================================================================================
63 /*
64  * Construct an empty InterfaceKinetics reaction mechanism.
65  * @param thermo This is an optional parameter that may be
66  * used to initialize the inherited Kinetics class with
67  * one ThermoPhase class object -> in other words it's
68  * useful for initialization of homogeneous kinetics
69  * mechanisms.
70  */
71 InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
72  Kinetics(),
73  m_redo_rates(false),
74  m_nirrev(0),
75  m_nrev(0),
76  m_surf(0),
77  m_integrator(0),
78  m_beta(0),
79  m_ctrxn(0),
80  m_ctrxn_ecdf(0),
81  m_StandardConc(0),
82  m_deltaG0(0),
83  m_ProdStanConcReac(0),
84  m_finalized(false),
85  m_has_coverage_dependence(false),
86  m_has_electrochem_rxns(false),
87  m_has_exchange_current_density_formulation(false),
88  m_phaseExistsCheck(false),
89  m_phaseExists(0),
90  m_phaseIsStable(0),
91  m_rxnPhaseIsReactant(0),
92  m_rxnPhaseIsProduct(0),
93  m_ioFlag(0)
94 {
95  if (thermo != 0) {
96  addPhase(*thermo);
97  }
99  m_kdata->m_temp = 0.0;
100 }
101 //====================================================================================================================
102 /*
103  * Destructor
104  */
106 {
107  delete m_kdata;
108  if (m_integrator) {
109  delete m_integrator;
110  }
111  for (size_t i = 0; i < m_ii; i++) {
112  delete [] m_rxnPhaseIsReactant[i];
113  delete [] m_rxnPhaseIsProduct[i];
114  }
115 }
116 //====================================================================================================================
117 // Copy Constructor for the %InterfaceKinetics object.
118 /*
119  * Currently, this is not fully implemented. If called it will
120  * throw an exception.
121  */
123  Kinetics(),
124  m_redo_rates(false),
125  m_nirrev(0),
126  m_nrev(0),
127  m_surf(0),
128  m_integrator(0),
129  m_beta(0),
130  m_ctrxn(0),
131  m_ctrxn_ecdf(0),
132  m_StandardConc(0),
133  m_deltaG0(0),
134  m_ProdStanConcReac(0),
135  m_finalized(false),
136  m_has_coverage_dependence(false),
137  m_has_electrochem_rxns(false),
138  m_has_exchange_current_density_formulation(false),
139  m_phaseExistsCheck(false),
140  m_phaseExists(0),
141  m_phaseIsStable(0),
142  m_rxnPhaseIsReactant(0),
143  m_rxnPhaseIsProduct(0),
144  m_ioFlag(0)
145 {
147  m_kdata->m_temp = 0.0;
148  /*
149  * Call the assignment operator
150  */
151  *this = operator=(right);
152 }
153 //====================================================================================================================
154 // Assignment operator
155 /*
156  * This is NOT a virtual function.
157  *
158  * @param right Reference to %Kinetics object to be copied into the
159  * current one.
160  */
163 {
164  /*
165  * Check for self assignment.
166  */
167  if (this == &right) {
168  return *this;
169  }
170 
171  for (size_t i = 0; i < m_ii; i++) {
172  delete [] m_rxnPhaseIsReactant[i];
173  delete [] m_rxnPhaseIsProduct[i];
174  }
175 
176  Kinetics::operator=(right);
177 
178  m_grt = right.m_grt;
179  m_revindex = right.m_revindex;
180  m_rates = right.m_rates;
181  m_redo_rates = right.m_redo_rates;
182  m_index = right.m_index;
183  m_irrev = right.m_irrev;
184  m_rxnstoich = right.m_rxnstoich;
185  m_nirrev = right.m_nirrev;
186  m_nrev = right.m_nrev;
187  m_rrxn = right.m_rrxn;
188  m_prxn = right.m_prxn;
189  m_rxneqn = right.m_rxneqn;
190  *m_kdata = *right.m_kdata;
191  m_conc = right.m_conc;
192  m_mu0 = right.m_mu0;
193  m_phi = right.m_phi;
194  m_pot = right.m_pot;
195  m_rwork = right.m_rwork;
196  m_E = right.m_E;
197  m_surf = right.m_surf; //DANGER - shallow copy
198  m_integrator = right.m_integrator; //DANGER - shallow copy
199  m_beta = right.m_beta;
200  m_ctrxn = right.m_ctrxn;
201  m_ctrxn_ecdf = right.m_ctrxn_ecdf;
202  m_StandardConc = right.m_StandardConc;
203  m_deltaG0 = right.m_deltaG0;
204  m_ProdStanConcReac = right.m_ProdStanConcReac;
205  m_finalized = right.m_finalized;
212 
213  m_rxnPhaseIsReactant.resize(m_ii, 0);
214  m_rxnPhaseIsProduct.resize(m_ii, 0);
215  size_t np = nPhases();
216  for (size_t i = 0; i < m_ii; i++) {
217  m_rxnPhaseIsReactant[i] = new bool[np];
218  m_rxnPhaseIsProduct[i] = new bool[np];
219  for (size_t p = 0; p < np; p++) {
220  m_rxnPhaseIsReactant[i][p] = right.m_rxnPhaseIsReactant[i][p];
221  m_rxnPhaseIsProduct[i][p] = right.m_rxnPhaseIsProduct[i][p];
222  }
223  }
224 
225  m_ioFlag = right.m_ioFlag;
226 
227  return *this;
228 }
229 //====================================================================================================================
230 // Return the ID of the kinetics object
232 {
233  return cInterfaceKinetics;
234 }
235 //====================================================================================================================
237 {
238  return cInterfaceKinetics;
239 }
240 //====================================================================================================================
241 // Duplication routine for objects which inherit from Kinetics
242 /*
243  * This virtual routine can be used to duplicate %Kinetics objects
244  * inherited from %Kinetics even if the application only has
245  * a pointer to %Kinetics to work with.
246  *
247  * These routines are basically wrappers around the derived copy
248  * constructor.
249  *
250  * @param tpVector Vector of shallow pointers to ThermoPhase objects. this is the
251  * m_thermo vector within this object
252  */
253 Kinetics* InterfaceKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
254 {
255  InterfaceKinetics* iK = new InterfaceKinetics(*this);
256  iK->assignShallowPointers(tpVector);
257  return dynamic_cast<Kinetics*>(iK);
258 }
259 //====================================================================================================================
260 // Set the electric potential in the nth phase
261 /*
262  * @param n phase Index in this kinetics object.
263  * @param V Electric potential (volts)
264  */
266 {
268  m_redo_rates = true;
269 }
270 //====================================================================================================================
271 // Update properties that depend on temperature
272 /*
273  * This is called to update all of the properties that depend on temperature
274  *
275  * Current objects that this function updates
276  * m_kdata->m_logtemp
277  * m_kdata->m_rfn
278  * m_rates.
279  * updateKc();
280  */
282 {
287  m_redo_rates = true;
288  }
289  doublereal T = thermo(surfacePhaseIndex()).temperature();
290  m_redo_rates = true;
291  if (T != m_kdata->m_temp || m_redo_rates) {
292  m_kdata->m_logtemp = log(T);
296  }
299  }
300  m_kdata->m_temp = T;
301  updateKc();
302  m_kdata->m_ROP_ok = false;
303  m_redo_rates = false;
304  }
305 }
306 //====================================================================================================================
308 {
309  for (size_t n = 0; n < nPhases(); n++) {
310  if (thermo(n).electricPotential() != m_phi[n]) {
311  m_phi[n] = thermo(n).electricPotential();
312  m_redo_rates = true;
313  }
314  }
315 }
316 //====================================================================================================================
317 
318 
319 /**
320  * Update properties that depend on concentrations. This method
321  * fills out the array of generalized concentrations by calling
322  * method getActivityConcentrations for each phase, which classes
323  * representing phases should overload to return the appropriate
324  * quantities.
325  */
327 {
328  for (size_t n = 0; n < nPhases(); n++) {
329  /*
330  * We call the getActivityConcentrations function of each
331  * ThermoPhase class that makes up this kinetics object to
332  * obtain the generalized concentrations for species within that
333  * class. This is collected in the vector m_conc. m_start[]
334  * are integer indices for that vector denoting the start of the
335  * species for each phase.
336  */
338  }
339  m_kdata->m_ROP_ok = false;
340 }
341 
342 
343 // Get the vector of activity concentrations used in the kinetics object
344 /*
345  * @param conc (output) Vector of activity concentrations. Length is
346  * equal to the number of species in the kinetics object
347  */
349 {
350  _update_rates_C();
351  copy(m_conc.begin(), m_conc.end(), conc);
352 }
353 
354 
355 /**
356  * Update the equilibrium constants in molar units for all
357  * reversible reactions. Irreversible reactions have their
358  * equilibrium constant set to zero.
359  */
361 {
362  vector_fp& m_rkc = m_kdata->m_rkcn;
363  fill(m_rkc.begin(), m_rkc.end(), 0.0);
364 
365  //static vector_fp mu(nTotalSpecies());
366  if (m_nrev > 0) {
367 
368  size_t nsp, ik = 0;
369  doublereal rt = GasConstant*thermo(0).temperature();
370  doublereal rrt = 1.0 / rt;
371  size_t np = nPhases();
372  for (size_t n = 0; n < np; n++) {
374  nsp = thermo(n).nSpecies();
375  for (size_t k = 0; k < nsp; k++) {
376  m_mu0[ik] -= rt * thermo(n).logStandardConc(k);
377  m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
378  ik++;
379  }
380  }
381 
382  // compute Delta mu^0 for all reversible reactions
384  DATA_PTR(m_rkc));
385 
386  for (size_t i = 0; i < m_nrev; i++) {
387  size_t irxn = m_revindex[i];
388  if (irxn == npos || irxn >= nReactions()) {
389  throw CanteraError("InterfaceKinetics",
390  "illegal value: irxn = "+int2str(irxn));
391  }
392  m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
393  }
394  for (size_t i = 0; i != m_nirrev; ++i) {
395  m_rkc[ m_irrev[i] ] = 0.0;
396  }
397  }
398 }
399 //====================================================================================================================
400 
401 
402 
403 void InterfaceKinetics::checkPartialEquil()
404 {
405  vector_fp dmu(nTotalSpecies(), 0.0);
406  vector_fp frop(std::max<size_t>(nReactions(), 1), 0.0);
407  vector_fp rrop(std::max<size_t>(nReactions(), 1), 0.0);
408  vector_fp netrop(std::max<size_t>(nReactions(), 1), 0.0);
409  vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
410  if (m_nrev > 0) {
411  doublereal rt = GasConstant*thermo(0).temperature();
412  cout << "T = " << thermo(0).temperature() << " " << rt << endl;
413  size_t nsp, ik=0;
414  //doublereal rt = GasConstant*thermo(0).temperature();
415  // doublereal rrt = 1.0/rt;
416  doublereal delta;
417  for (size_t n = 0; n < nPhases(); n++) {
419  nsp = thermo(n).nSpecies();
420  for (size_t k = 0; k < nsp; k++) {
421  delta = Faraday * m_phi[n] * thermo(n).charge(k);
422  //cout << thermo(n).speciesName(k) << " " << (delta+dmu[ik])/rt << " " << dmu[ik]/rt << endl;
423  dmu[ik] += delta;
424  ik++;
425  }
426  }
427 
428  // compute Delta mu^ for all reversible reactions
433  for (size_t i = 0; i < m_nrev; i++) {
434  size_t irxn = m_revindex[i];
435  cout << "Reaction " << reactionString(irxn)
436  << " " << rmu[irxn]/rt << endl;
437  printf("%12.6e %12.6e %12.6e %12.6e \n",
438  frop[irxn], rrop[irxn], netrop[irxn],
439  netrop[irxn]/(frop[irxn] + rrop[irxn]));
440  }
441  }
442 }
443 
444 
445 /**
446  * Get the equilibrium constants of all reactions, whether
447  * reversible or not.
448  */
450 {
451  size_t ik=0;
452  doublereal rt = GasConstant*thermo(0).temperature();
453  doublereal rrt = 1.0/rt;
454  for (size_t n = 0; n < nPhases(); n++) {
456  size_t nsp = thermo(n).nSpecies();
457  for (size_t k = 0; k < nsp; k++) {
458  m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
459  m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
460  ik++;
461  }
462  }
463 
464  fill(kc, kc + m_ii, 0.0);
465 
467 
468  for (size_t i = 0; i < m_ii; i++) {
469  kc[i] = exp(-kc[i]*rrt);
470  }
471 }
472 
473 void InterfaceKinetics::getExchangeCurrentQuantities()
474 {
475  /*
476  * First collect vectors of the standard Gibbs free energies of the
477  * species and the standard concentrations
478  * - m_mu0
479  * - m_logStandardConc
480  */
481  size_t ik = 0;
482 
483  for (size_t n = 0; n < nPhases(); n++) {
485  size_t nsp = thermo(n).nSpecies();
486  for (size_t k = 0; k < nsp; k++) {
487  m_StandardConc[ik] = thermo(n).standardConcentration(k);
488  ik++;
489  }
490  }
491 
493 
494  for (size_t i = 0; i < m_ii; i++) {
495  m_ProdStanConcReac[i] = 1.0;
496  }
497 
498  m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac));
499 
500 }
501 
502 // Returns the Species creation rates [kmol/m^2/s].
503 /*
504  * Return the species
505  * creation rates in array cdot, which must be
506  * dimensioned at least as large as the total number of
507  * species in all phases of the kinetics
508  * model
509  *
510  * @param cdot Vector containing the creation rates.
511  * length = m_kk. units = kmol/m^2/s
512  */
514 {
515  updateROP();
517  &m_kdata->m_ropr[0], cdot);
518 }
519 
520 // Return the Species destruction rates [kmol/m^2/s].
521 /*
522  * Return the species destruction rates in array ddot, which must be
523  * dimensioned at least as large as the total number of
524  * species in all phases of the kinetics model
525  */
527 {
528  updateROP();
530  &m_kdata->m_ropr[0], ddot);
531 }
532 
533 // Return the species net production rates
534 /*
535  * Species net production rates [kmol/m^2/s]. Return the species
536  * net production rates (creation - destruction) in array
537  * wdot, which must be dimensioned at least as large as the
538  * total number of species in all phases of the kinetics
539  * model
540  *
541  * @param net Vector of species production rates.
542  * units kmol m-d s-1, where d is dimension.
543  */
545 {
546  updateROP();
548  &m_kdata->m_ropnet[0],
549  net);
550 }
551 
552 //====================================================================================================================
553 // Apply corrections for interfacial charge transfer reactions
554 /*
555  * For reactions that transfer charge across a potential difference,
556  * the activation energies are modified by the potential difference.
557  * (see, for example, ...). This method applies this correction.
558  *
559  * @param kf Vector of forward reaction rate constants on which to have
560  * the correction applied
561  */
563 {
564  // compute the electrical potential energy of each species
565  size_t ik = 0;
566  for (size_t n = 0; n < nPhases(); n++) {
567  size_t nsp = thermo(n).nSpecies();
568  for (size_t k = 0; k < nsp; k++) {
569  m_pot[ik] = Faraday*thermo(n).charge(k)*m_phi[n];
570  ik++;
571  }
572  }
573 
574  // Compute the change in electrical potential energy for each
575  // reaction. This will only be non-zero if a potential
576  // difference is present.
578 
579  // Modify the reaction rates. Only modify those with a
580  // non-zero activation energy. Below we decrease the
581  // activation energy below zero but in some debug modes
582  // we print out a warning message about this.
583  /*
584  * NOTE, there is some discussion about this point.
585  * Should we decrease the activation energy below zero?
586  * I don't think this has been decided in any definitive way.
587  * The treatment below is numerically more stable, however.
588  */
589  doublereal eamod;
590 #ifdef DEBUG_KIN_MODE
591  doublereal ea;
592 #endif
593  for (size_t i = 0; i < m_beta.size(); i++) {
594  size_t irxn = m_ctrxn[i];
595  eamod = m_beta[i]*m_rwork[irxn];
596  // if (eamod != 0.0 && m_E[irxn] != 0.0) {
597  if (eamod != 0.0) {
598 #ifdef DEBUG_KIN_MODE
599  ea = GasConstant * m_E[irxn];
600  if (eamod + ea < 0.0) {
601  writelog("Warning: act energy mod too large!\n");
602  writelog(" Delta phi = "+fp2str(m_rwork[irxn]/Faraday)+"\n");
603  writelog(" Delta Ea = "+fp2str(eamod)+"\n");
604  writelog(" Ea = "+fp2str(ea)+"\n");
605  for (n = 0; n < np; n++) {
606  writelog("Phase "+int2str(n)+": phi = "
607  +fp2str(m_phi[n])+"\n");
608  }
609  }
610 #endif
611  doublereal rt = GasConstant*thermo(0).temperature();
612  doublereal rrt = 1.0/rt;
613  kf[irxn] *= exp(-eamod*rrt);
614  }
615  }
616 }
617 //====================================================================================================================
619 {
620  getExchangeCurrentQuantities();
621  doublereal rt = GasConstant*thermo(0).temperature();
622  doublereal rrt = 1.0/rt;
623  for (size_t i = 0; i < m_ctrxn.size(); i++) {
624  size_t irxn = m_ctrxn[i];
625  int iECDFormulation = m_ctrxn_ecdf[i];
626  if (iECDFormulation) {
627  double tmp = exp(- m_beta[i] * m_deltaG0[irxn] * rrt);
628  double tmp2 = m_ProdStanConcReac[irxn];
629  tmp *= 1.0 / tmp2 / Faraday;
630  kfwd[irxn] *= tmp;
631  }
632  }
633 
634 
635 }
636 //====================================================================================================================
637 /**
638  * Update the rates of progress of the reactions in the reaction
639  * mechanism. This routine operates on internal data.
640  */
642 {
643 
644  updateROP();
645 
646  const vector_fp& rf = m_kdata->m_rfn;
647 
648  // copy rate coefficients into kfwd
649  copy(rf.begin(), rf.end(), kfwd);
650 
651  // multiply by perturbation factor
652  multiply_each(kfwd, kfwd + nReactions(), m_perturb.begin());
653 
654 }
655 //====================================================================================================================
656 
657 /**
658  * Update the rates of progress of the reactions in the reaction
659  * mechanism. This routine operates on internal data.
660  */
661 void InterfaceKinetics::getRevRateConstants(doublereal* krev, bool doIrreversible)
662 {
663  getFwdRateConstants(krev);
664  if (doIrreversible) {
665  doublereal* tmpKc = DATA_PTR(m_kdata->m_ropnet);
667  for (size_t i = 0; i < m_ii; i++) {
668  krev[i] /= tmpKc[i];
669  }
670  } else {
671  const vector_fp& rkc = m_kdata->m_rkcn;
672  multiply_each(krev, krev + nReactions(), rkc.begin());
673  }
674 }
675 //====================================================================================================================
676 
678 {
679  copy(m_E.begin(), m_E.end(), E);
680 }
681 //====================================================================================================================
682 /**
683  * Update the rates of progress of the reactions in the reaction
684  * mechanism. This routine operates on internal data.
685  */
687 {
688 
689  _update_rates_T();
690  _update_rates_C();
691 
692  if (m_kdata->m_ROP_ok) {
693  return;
694  }
695 
696  const vector_fp& rf = m_kdata->m_rfn;
697  const vector_fp& m_rkc = m_kdata->m_rkcn;
698  vector_fp& ropf = m_kdata->m_ropf;
699  vector_fp& ropr = m_kdata->m_ropr;
700  vector_fp& ropnet = m_kdata->m_ropnet;
701 
702  // copy rate coefficients into ropf
703  copy(rf.begin(), rf.end(), ropf.begin());
704 
705  // multiply by perturbation factor
706  multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
707 
708  // copy the forward rates to the reverse rates
709  copy(ropf.begin(), ropf.end(), ropr.begin());
710 
711  // for reverse rates computed from thermochemistry, multiply
712  // the forward rates copied into m_ropr by the reciprocals of
713  // the equilibrium constants
714  multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());
715 
716  // multiply ropf by concentration products
718  //m_reactantStoich.multiply(m_conc.begin(), ropf.begin());
719 
720  // for reversible reactions, multiply ropr by concentration
721  // products
723  DATA_PTR(ropr));
724  //m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
725 
726  // do global reactions
727  //m_globalReactantStoich.power(m_conc.begin(), ropf.begin());
728 
729  for (size_t j = 0; j != m_ii; ++j) {
730  ropnet[j] = ropf[j] - ropr[j];
731  }
732 
733 
734  /*
735  * For reactions involving multiple phases, we must check that the phase
736  * being consumed actually exists. This is particularly important for
737  * phases that are stoichiometric phases containing one species with a unity activity
738  */
739  if (m_phaseExistsCheck) {
740  for (size_t j = 0; j != m_ii; ++j) {
741  if ((ropr[j] > ropf[j]) && (ropr[j] > 0.0)) {
742  for (size_t p = 0; p < nPhases(); p++) {
743  if (m_rxnPhaseIsProduct[j][p]) {
744  if (! m_phaseExists[p]) {
745  ropnet[j] = 0.0;
746  ropr[j] = ropf[j];
747  if (ropf[j] > 0.0) {
748  for (size_t rp = 0; rp < nPhases(); rp++) {
749  if (m_rxnPhaseIsReactant[j][rp]) {
750  if (! m_phaseExists[rp]) {
751  ropnet[j] = 0.0;
752  ropr[j] = ropf[j] = 0.0;
753  }
754  }
755  }
756  }
757  }
758  }
759  if (m_rxnPhaseIsReactant[j][p]) {
760  if (! m_phaseIsStable[p]) {
761  ropnet[j] = 0.0;
762  ropr[j] = ropf[j];
763  }
764  }
765  }
766  } else if ((ropf[j] > ropr[j]) && (ropf[j] > 0.0)) {
767  for (size_t p = 0; p < nPhases(); p++) {
768  if (m_rxnPhaseIsReactant[j][p]) {
769  if (! m_phaseExists[p]) {
770  ropnet[j] = 0.0;
771  ropf[j] = ropr[j];
772  if (ropf[j] > 0.0) {
773  for (size_t rp = 0; rp < nPhases(); rp++) {
774  if (m_rxnPhaseIsProduct[j][rp]) {
775  if (! m_phaseExists[rp]) {
776  ropnet[j] = 0.0;
777  ropf[j] = ropr[j] = 0.0;
778  }
779  }
780  }
781  }
782  }
783  }
784  if (m_rxnPhaseIsProduct[j][p]) {
785  if (! m_phaseIsStable[p]) {
786  ropnet[j] = 0.0;
787  ropf[j] = ropr[j];
788  }
789  }
790  }
791  }
792  }
793  }
794 
795  m_kdata->m_ROP_ok = true;
796 }
797 
798 #ifdef KINETICS_WITH_INTERMEDIATE_ZEROED_PHASES
799 //=================================================================================================
800 InterfaceKinetics::adjustRatesForIntermediatePhases()
801 {
802  doublereal sFac = 1.0;
803 
804  vector_fp& ropf = m_kdata->m_ropf;
805  vector_fp& ropr = m_kdata->m_ropr;
806  vector_fp& ropnet = m_kdata->m_ropnet;
807 
808  getCreatingRates(DATA_PTR(m_speciestmpP));
809  getDestructionRates(DATA_PTR(m_speciestmpD));
810 
811  for (iphase = 0; iphase < nphases; iphase++) {
812  if (m_intermediatePhases(iphase)) {
813  for (isp = 0; isp < nspecies; isp++) {
814  if (m_speciesTmpD[ispI] > m_speciesTmpP[I]) {
815  sFac = m_speciesTmpD[ispI]/ m_speciesTmpP[I];
816  }
817  // Loop over reactions that are reactants for the species in the phase
818  // reducing their rates.
819 
820 
821  }
822  }
823 
824  }
825 
826 }
827 #endif
828 //=================================================================================================
829 //=================================================================================================
830 /*
831  *
832  * getDeltaGibbs():
833  *
834  * Return the vector of values for the reaction gibbs free energy
835  * change
836  * These values depend upon the concentration
837  * of the ideal gas.
838  *
839  * units = J kmol-1
840  */
841 void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG)
842 {
843  /*
844  * Get the chemical potentials of the species in the
845  * ideal gas solution.
846  */
847  for (size_t n = 0; n < nPhases(); n++) {
849  }
850  //for (n = 0; n < m_grt.size(); n++) {
851  // cout << n << "G_RT = " << m_grt[n] << endl;
852  //}
853 
854  /*
855  * Use the stoichiometric manager to find deltaG for each
856  * reaction.
857  */
858  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
859 }
860 //=================================================================================================
861 // Return the vector of values for the reaction electrochemical free energy change.
862 /*
863  * These values depend upon the concentration of the solution and
864  * the voltage of the phases
865  *
866  * units = J kmol-1
867  *
868  * @param deltaM Output vector of deltaM's for reactions
869  * Length: m_ii.
870  */
872 {
873  /*
874  * Get the chemical potentials of the species in the
875  * ideal gas solution.
876  */
877  size_t np = nPhases();
878  for (size_t n = 0; n < np; n++) {
880  }
881  /*
882  * Use the stoichiometric manager to find deltaG for each
883  * reaction.
884  */
885  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaM);
886 }
887 //=================================================================================================
888 /*
889  *
890  * getDeltaEnthalpy():
891  *
892  * Return the vector of values for the reactions change in
893  * enthalpy.
894  * These values depend upon the concentration
895  * of the solution.
896  *
897  * units = J kmol-1
898  */
899 void InterfaceKinetics::getDeltaEnthalpy(doublereal* deltaH)
900 {
901  /*
902  * Get the partial molar enthalpy of all species in the
903  * ideal gas.
904  */
905  for (size_t n = 0; n < nPhases(); n++) {
907  }
908  /*
909  * Use the stoichiometric manager to find deltaG for each
910  * reaction.
911  */
912  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
913 }
914 
915 
916 // Return the vector of values for the change in
917 // entropy due to each reaction
918 /*
919  * These values depend upon the concentration
920  * of the solution.
921  *
922  * units = J kmol-1 Kelvin-1
923  *
924  * @param deltaS vector of Enthalpy changes
925  * Length = m_ii, number of reactions
926  *
927  */
928 void InterfaceKinetics::getDeltaEntropy(doublereal* deltaS)
929 {
930  /*
931  * Get the partial molar entropy of all species in all of
932  * the phases
933  */
934  for (size_t n = 0; n < nPhases(); n++) {
936  }
937  /*
938  * Use the stoichiometric manager to find deltaS for each
939  * reaction.
940  */
941  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
942 }
943 
944 /**
945  *
946  * getDeltaSSGibbs():
947  *
948  * Return the vector of values for the reaction
949  * standard state gibbs free energy change.
950  * These values don't depend upon the concentration
951  * of the solution.
952  *
953  * units = J kmol-1
954  */
955 void InterfaceKinetics::getDeltaSSGibbs(doublereal* deltaG)
956 {
957  /*
958  * Get the standard state chemical potentials of the species.
959  * This is the array of chemical potentials at unit activity
960  * We define these here as the chemical potentials of the pure
961  * species at the temperature and pressure of the solution.
962  */
963  for (size_t n = 0; n < nPhases(); n++) {
965  }
966  /*
967  * Use the stoichiometric manager to find deltaG for each
968  * reaction.
969  */
970  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaG);
971 }
972 
973 /**
974  *
975  * getDeltaSSEnthalpy():
976  *
977  * Return the vector of values for the change in the
978  * standard state enthalpies of reaction.
979  * These values don't depend upon the concentration
980  * of the solution.
981  *
982  * units = J kmol-1
983  */
985 {
986  /*
987  * Get the standard state enthalpies of the species.
988  * This is the array of chemical potentials at unit activity
989  * We define these here as the enthalpies of the pure
990  * species at the temperature and pressure of the solution.
991  */
992  for (size_t n = 0; n < nPhases(); n++) {
994  }
995  doublereal RT = thermo().temperature() * GasConstant;
996  for (size_t k = 0; k < m_kk; k++) {
997  m_grt[k] *= RT;
998  }
999  /*
1000  * Use the stoichiometric manager to find deltaG for each
1001  * reaction.
1002  */
1003  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaH);
1004 }
1005 
1006 /*********************************************************************
1007  *
1008  * getDeltaSSEntropy():
1009  *
1010  * Return the vector of values for the change in the
1011  * standard state entropies for each reaction.
1012  * These values don't depend upon the concentration
1013  * of the solution.
1014  *
1015  * units = J kmol-1 Kelvin-1
1016  */
1018 {
1019  /*
1020  * Get the standard state entropy of the species.
1021  * We define these here as the entropies of the pure
1022  * species at the temperature and pressure of the solution.
1023  */
1024  for (size_t n = 0; n < nPhases(); n++) {
1026  }
1027  doublereal R = GasConstant;
1028  for (size_t k = 0; k < m_kk; k++) {
1029  m_grt[k] *= R;
1030  }
1031  /*
1032  * Use the stoichiometric manager to find deltaS for each
1033  * reaction.
1034  */
1035  m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_grt), deltaS);
1036 }
1037 
1038 //====================================================================================================================
1039 /**
1040  * Add a single reaction to the mechanism. This routine
1041  * must be called after init() and before finalize().
1042  * This function branches on the types of reactions allowed
1043  * by the interfaceKinetics manager in order to install
1044  * the reaction correctly in the manager.
1045  * The manager allows the following reaction types
1046  * Elementary
1047  * Surface
1048  * Global
1049  * There is no difference between elementary and surface
1050  * reactions.
1051  */
1052 void InterfaceKinetics::addReaction(ReactionData& r)
1053 {
1054 
1055  /*
1056  * Install the rate coefficient for the current reaction
1057  * in the appropriate data structure.
1058  */
1059  addElementaryReaction(r);
1060  /*
1061  * Add the reactants and products for m_ropnet;the current reaction
1062  * to the various stoichiometric coefficient arrays.
1063  */
1064  installReagents(r);
1065  /*
1066  * Save the reaction and product groups, which are
1067  * part of the ReactionData class, in this class.
1068  * They aren't used for anything but reaction path
1069  * analysis.
1070  */
1071  //installGroups(reactionNumber(), r.rgroups, r.pgroups);
1072  /*
1073  * Increase the internal number of reactions, m_ii, by one.
1074  * increase the size of m_perturb by one as well.
1075  */
1077  m_rxneqn.push_back(r.equation);
1078 
1079  m_rxnPhaseIsReactant.resize(m_ii, 0);
1080  m_rxnPhaseIsProduct.resize(m_ii, 0);
1081 
1082  size_t np = nPhases();
1083  size_t i = m_ii - 1;
1084  m_rxnPhaseIsReactant[i] = new bool[np];
1085  m_rxnPhaseIsProduct[i] = new bool[np];
1086 
1087  for (size_t p = 0; p < np; p++) {
1088  m_rxnPhaseIsReactant[i][p] = false;
1089  m_rxnPhaseIsProduct[i][p] = false;
1090  }
1091 
1092  const std::vector<size_t>& vr = reactants(i);
1093  for (size_t ik = 0; ik < vr.size(); ik++) {
1094  size_t k = vr[ik];
1095  size_t p = speciesPhaseIndex(k);
1096  m_rxnPhaseIsReactant[i][p] = true;
1097  }
1098  const std::vector<size_t>& vp = products(i);
1099  for (size_t ik = 0; ik < vp.size(); ik++) {
1100  size_t k = vp[ik];
1101  size_t p = speciesPhaseIndex(k);
1102  m_rxnPhaseIsProduct[i][p] = true;
1103  }
1104 }
1105 //====================================================================================================================
1106 void InterfaceKinetics::addElementaryReaction(ReactionData& r)
1107 {
1108  // install rate coeff calculator
1109  vector_fp& rp = r.rateCoeffParameters;
1110  size_t ncov = r.cov.size();
1111  if (ncov > 3) {
1113  }
1114  for (size_t m = 0; m < ncov; m++) {
1115  rp.push_back(r.cov[m]);
1116  }
1117 
1118  /*
1119  * Temporarily change the reaction rate coefficient type to surface arrhenius.
1120  * This is what is expected. We'll handle exchange current types below by hand.
1121  */
1122  int reactionRateCoeffType_orig = r.rateCoeffType;
1123  if (r.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
1124  r.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
1125  }
1126  if (r.rateCoeffType == ARRHENIUS_REACTION_RATECOEFF_TYPE) {
1127  r.rateCoeffType = SURF_ARRHENIUS_REACTION_RATECOEFF_TYPE;
1128  }
1129  /*
1130  * Install the reaction rate into the vector of reactions handled by this class
1131  */
1132  size_t iloc = m_rates.install(reactionNumber(), r);
1133 
1134  /*
1135  * Change the reaction rate coefficient type back to its original value
1136  */
1137  r.rateCoeffType = reactionRateCoeffType_orig;
1138 
1139  // store activation energy
1140  m_E.push_back(r.rateCoeffParameters[2]);
1141 
1142  if (r.beta > 0.0) {
1143  m_has_electrochem_rxns = true;
1144  m_beta.push_back(r.beta);
1145  m_ctrxn.push_back(reactionNumber());
1146  if (r.rateCoeffType == EXCHANGE_CURRENT_REACTION_RATECOEFF_TYPE) {
1148  m_ctrxn_ecdf.push_back(1);
1149  } else {
1150  m_ctrxn_ecdf.push_back(0);
1151  }
1152  }
1153 
1154  // add constant term to rate coeff value vector
1155  m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
1156  registerReaction(reactionNumber(), ELEMENTARY_RXN, iloc);
1157 }
1158 //====================================================================================================================
1159 
1160 void InterfaceKinetics::setIOFlag(int ioFlag)
1161 {
1162  m_ioFlag = ioFlag;
1163  if (m_integrator) {
1164  m_integrator->setIOFlag(ioFlag);
1165  }
1166 }
1167 
1168 // void InterfaceKinetics::
1169 // addGlobalReaction(const ReactionData& r) {
1170 
1171 // int iloc;
1172 // // install rate coeff calculator
1173 // vector_fp rp = r.rateCoeffParameters;
1174 // int ncov = r.cov.size();
1175 // for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]);
1176 // iloc = m_rates.install( reactionNumber(),
1177 // r.rateCoeffType, rp.size(),
1178 // rp.begin() );
1179 // // store activation energy
1180 // m_E.push_back(r.rateCoeffParameters[2]);
1181 // // add constant term to rate coeff value vector
1182 // m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);
1183 
1184 // int nr = r.order.size();
1185 // vector_fp ordr(nr);
1186 // for (int n = 0; n < nr; n++) {
1187 // ordr[n] = r.order[n] - r.rstoich[n];
1188 // }
1189 // m_globalReactantStoich.add( reactionNumber(),
1190 // r.reactants, ordr);
1191 
1192 // registerReaction( reactionNumber(), GLOBAL_RXN, iloc);
1193 // }
1194 
1195 
1196 void InterfaceKinetics::installReagents(const ReactionData& r)
1197 {
1198 
1199  size_t n, ns, m;
1200  doublereal nsFlt;
1201  /*
1202  * extend temporary storage by one for this rxn.
1203  */
1204  m_kdata->m_ropf.push_back(0.0);
1205  m_kdata->m_ropr.push_back(0.0);
1206  m_kdata->m_ropnet.push_back(0.0);
1207  m_kdata->m_rkcn.push_back(0.0);
1208 
1209  /*
1210  * Obtain the current reaction index for the reaction that we
1211  * are adding. The first reaction is labeled 0.
1212  */
1213  size_t rnum = reactionNumber();
1214 
1215  // vectors rk and pk are lists of species numbers, with
1216  // repeated entries for species with stoichiometric
1217  // coefficients > 1. This allows the reaction to be defined
1218  // with unity reaction order for each reactant, and so the
1219  // faster method 'multiply' can be used to compute the rate of
1220  // progress instead of 'power'.
1221 
1222  std::vector<size_t> rk;
1223  size_t nr = r.reactants.size();
1224  for (n = 0; n < nr; n++) {
1225  nsFlt = r.rstoich[n];
1226  ns = (size_t) nsFlt;
1227  if ((doublereal) ns != nsFlt) {
1228  if (ns < 1) {
1229  ns = 1;
1230  }
1231  }
1232  /*
1233  * Add to m_rrxn. m_rrxn is a vector of maps. m_rrxn has a length
1234  * equal to the total number of species for each species, there
1235  * exists a map, with the reaction number being the key, and the
1236  * reactant stoichiometric coefficient being the value.
1237  */
1238  m_rrxn[r.reactants[n]][rnum] = nsFlt;
1239  for (m = 0; m < ns; m++) {
1240  rk.push_back(r.reactants[n]);
1241  }
1242  }
1243  /*
1244  * Now that we have rk[], we add it into the vector<vector_int> m_reactants
1245  * in the rnum index spot. Thus m_reactants[rnum] yields a vector
1246  * of reactants for the rnum'th reaction
1247  */
1248  m_reactants.push_back(rk);
1249  std::vector<size_t> pk;
1250  size_t np = r.products.size();
1251  for (n = 0; n < np; n++) {
1252  nsFlt = r.pstoich[n];
1253  ns = (size_t) nsFlt;
1254  if ((doublereal) ns != nsFlt) {
1255  if (ns < 1) {
1256  ns = 1;
1257  }
1258  }
1259  /*
1260  * Add to m_prxn. m_prxn is a vector of maps. m_prxn has a length
1261  * equal to the total number of species for each species, there
1262  * exists a map, with the reaction number being the key, and the
1263  * product stoichiometric coefficient being the value.
1264  */
1265  m_prxn[r.products[n]][rnum] = nsFlt;
1266  for (m = 0; m < ns; m++) {
1267  pk.push_back(r.products[n]);
1268  }
1269  }
1270  /*
1271  * Now that we have pk[], we add it into the vector<vector_int> m_products
1272  * in the rnum index spot. Thus m_products[rnum] yields a vector
1273  * of products for the rnum'th reaction
1274  */
1275  m_products.push_back(pk);
1276  /*
1277  * Add this reaction to the stoichiometric coefficient manager. This
1278  * calculates rates of species production from reaction rates of
1279  * progress.
1280  */
1281  m_rxnstoich.add(reactionNumber(), r);
1282  /*
1283  * register reaction in lists of reversible and irreversible rxns.
1284  */
1285  if (r.reversible) {
1286  m_revindex.push_back(reactionNumber());
1287  m_nrev++;
1288  } else {
1289  m_irrev.push_back(reactionNumber());
1290  m_nirrev++;
1291  }
1292 }
1293 //===============================================================================================
1295 {
1296  Kinetics::addPhase(thermo);
1297  m_phaseExists.push_back(true);
1298  m_phaseIsStable.push_back(true);
1299 }
1300 //================================================================================================
1301 /**
1302  * Prepare the class for the addition of reactions. This function
1303  * must be called after instantiation of the class, but before
1304  * any reactions are actually added to the mechanism.
1305  * This function calculates m_kk the number of species in all
1306  * phases participating in the reaction mechanism. We don't know
1307  * m_kk previously, before all phases have been added.
1308  */
1310 {
1311  m_kk = 0;
1312  for (size_t n = 0; n < nPhases(); n++) {
1313  m_kk += thermo(n).nSpecies();
1314  }
1315  m_rrxn.resize(m_kk);
1316  m_prxn.resize(m_kk);
1317  m_conc.resize(m_kk);
1318  m_mu0.resize(m_kk);
1319  m_grt.resize(m_kk);
1320  m_pot.resize(m_kk, 0.0);
1321  m_phi.resize(nPhases(), 0.0);
1322 }
1323 //================================================================================================
1324 /**
1325  * Finish adding reactions and prepare for use. This function
1326  * must be called after all reactions are entered into the mechanism
1327  * and before the mechanism is used to calculate reaction rates.
1328  *
1329  * Here, we resize work arrays based on the number of reactions,
1330  * since we don't know this number up to now.
1331  */
1333 {
1335  size_t safe_reaction_size = std::max<size_t>(nReactions(), 1);
1336  m_rwork.resize(safe_reaction_size);
1337  size_t ks = reactionPhaseIndex();
1338  if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
1339  "no surface phase is present.");
1340  m_surf = (SurfPhase*)&thermo(ks);
1341  if (m_surf->nDim() != 2)
1342  throw CanteraError("InterfaceKinetics::finalize",
1343  "expected interface dimension = 2, but got dimension = "
1344  +int2str(m_surf->nDim()));
1345 
1346  m_StandardConc.resize(m_kk, 0.0);
1347  m_deltaG0.resize(safe_reaction_size, 0.0);
1348  m_ProdStanConcReac.resize(safe_reaction_size, 0.0);
1349 
1350  if (m_thermo.size() != m_phaseExists.size()) {
1351  throw CanteraError("InterfaceKinetics::finalize", "internal error");
1352  }
1353 
1354  // Guarantee that these arrays can be converted to double* even in the
1355  // special case where there are no reactions defined.
1356  if (!m_ii) {
1357  m_perturb.resize(1, 1.0);
1358  }
1359 
1360  m_finalized = true;
1361 }
1362 
1363 doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
1364 {
1365  for (size_t i = 0; i < m_ctrxn.size(); i++) {
1366  if (m_ctrxn[i] == irxn) {
1367  return m_beta[i];
1368  }
1369  }
1370  return 0.0;
1371 }
1372 
1373 //================================================================================================
1375 {
1376  return (m_finalized);
1377 }
1378 //================================================================================================
1379 // Advance the surface coverages in time
1380 /*
1381  * @param tstep Time value to advance the surface coverages
1382  */
1384 advanceCoverages(doublereal tstep)
1385 {
1386  if (m_integrator == 0) {
1387  vector<InterfaceKinetics*> k;
1388  k.push_back(this);
1389  m_integrator = new ImplicitSurfChem(k);
1391  }
1392  m_integrator->integrate(0.0, tstep);
1393  delete m_integrator;
1394  m_integrator = 0;
1395 }
1396 //================================================================================================
1397 // Solve for the pseudo steady-state of the surface problem
1398 /*
1399  * Solve for the steady state of the surface problem.
1400  * This is the same thing as the advanceCoverages() function,
1401  * but at infinite times.
1402  *
1403  * Note, a direct solve is carried out under the hood here,
1404  * to reduce the computational time.
1405  *
1406  * the integrator object is saved between calls to
1407  * reduce the computational cost of repeated calls.
1408  */
1410 solvePseudoSteadyStateProblem(int ifuncOverride, doublereal timeScaleOverride)
1411 {
1412  // create our own solver object
1413  if (m_integrator == 0) {
1414  vector<InterfaceKinetics*> k;
1415  k.push_back(this);
1416  m_integrator = new ImplicitSurfChem(k);
1418  }
1419  m_integrator->setIOFlag(m_ioFlag);
1420  /*
1421  * New direct method to go here
1422  */
1423  m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride);
1424 }
1425 //================================================================================================
1426 
1427 void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists)
1428 {
1429  if (iphase >= m_thermo.size()) {
1430  throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds");
1431  }
1432  if (exists) {
1433  if (!m_phaseExists[iphase]) {
1435  m_phaseExists[iphase] = true;
1436  }
1437  m_phaseIsStable[iphase] = true;
1438  } else {
1439  if (m_phaseExists[iphase]) {
1441  m_phaseExists[iphase] = false;
1442  }
1443  m_phaseIsStable[iphase] = false;
1444  }
1445 }
1446 //================================================================================================
1447 // Gets the phase existence int for the ith phase
1448 /*
1449  * @param iphase Phase Id
1450  *
1451  * @return Returns the int specifying whether the kinetics object thinks the phase exists
1452  * or not. If it exists, then species in that phase can be a reactant in reactions.
1453  */
1454 int InterfaceKinetics::phaseExistence(const int iphase) const
1455 {
1456  if (iphase < 0 || iphase >= (int) m_thermo.size()) {
1457  throw CanteraError("InterfaceKinetics:phaseExistence()", "out of bounds");
1458  }
1459  return m_phaseExists[iphase];
1460 }
1461 //================================================================================================
1462 // Gets the phase stability int for the ith phase
1463 /*
1464  * @param iphase Phase Id
1465  *
1466  * @return Returns the int specifying whether the kinetics object thinks the phase is stable
1467  * with nonzero mole numbers.
1468  * If it stable, then the kinetics object will allow for rates of production of
1469  * of species in that phase that are positive.
1470  */
1471 int InterfaceKinetics::phaseStability(const int iphase) const
1472 {
1473  if (iphase < 0 || iphase >= (int) m_thermo.size()) {
1474  throw CanteraError("InterfaceKinetics:phaseStability()", "out of bounds");
1475  }
1476  return m_phaseIsStable[iphase];
1477 }
1478 //================================================================================================
1479 
1480 void InterfaceKinetics::setPhaseStability(const int iphase, const int isStable)
1481 {
1482  if (iphase < 0 || iphase >= (int) m_thermo.size()) {
1483  throw CanteraError("InterfaceKinetics:setPhaseStability", "out of bounds");
1484  }
1485  if (isStable) {
1486  m_phaseIsStable[iphase] = true;
1487  } else {
1488  m_phaseIsStable[iphase] = false;
1489  }
1490 }
1491 
1492 //================================================================================================
1494 {
1495  m_rwork.resize(std::max<size_t>(nReactions(), 1));
1496  size_t ks = reactionPhaseIndex();
1497  if (ks == npos) throw CanteraError("EdgeKinetics::finalize",
1498  "no edge phase is present.");
1499  m_surf = (SurfPhase*)&thermo(ks);
1500  if (m_surf->nDim() != 1)
1501  throw CanteraError("EdgeKinetics::finalize",
1502  "expected interface dimension = 1, but got dimension = "
1503  +int2str(m_surf->nDim()));
1504 
1505  // Guarantee that these arrays can be converted to double* even in the
1506  // special case where there are no reactions defined.
1507  if (!m_ii) {
1508  m_perturb.resize(1, 1.0);
1509  }
1510 
1511  m_finalized = true;
1512 }
1513 //================================================================================================
1514 }
1515 
1516