Cantera  2.4.0
ImplicitSurfChem.cpp
Go to the documentation of this file.
1 /**
2  * @file ImplicitSurfChem.cpp
3  * Definitions for the implicit integration of surface site density equations
4  * (see \ref kineticsmgr and class
5  * \link Cantera::ImplicitSurfChem ImplicitSurfChem\endlink).
6  */
7 
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at http://www.cantera.org/license.txt for license and copyright information.
10 
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
21  m_nv(0),
22  m_numTotalBulkSpecies(0),
23  m_numTotalSpecies(0),
24  m_atol(1.e-14),
25  m_rtol(1.e-7),
26  m_maxstep(0.0),
27  m_mediumSpeciesStart(-1),
28  m_bulkSpeciesStart(-1),
29  m_surfSpeciesStart(-1),
30  m_commonTempPressForPhases(true),
31  m_ioFlag(0)
32 {
33  size_t ntmax = 0;
34  size_t kinSpIndex = 0;
35  // Loop over the number of surface kinetics objects
36  for (size_t n = 0; n < k.size(); n++) {
37  InterfaceKinetics* kinPtr = k[n];
38  m_vecKinPtrs.push_back(kinPtr);
39  size_t ns = k[n]->surfacePhaseIndex();
40  if (ns == npos) {
41  throw CanteraError("ImplicitSurfChem",
42  "kinetics manager contains no surface phase");
43  }
44  m_surfindex.push_back(ns);
45  m_surf.push_back((SurfPhase*)&k[n]->thermo(ns));
46  size_t nsp = m_surf.back()->nSpecies();
47  m_nsp.push_back(nsp);
48  m_nv += m_nsp.back();
49  size_t nt = k[n]->nTotalSpecies();
50  ntmax = std::max(nt, ntmax);
51  m_specStartIndex.push_back(kinSpIndex);
52  kinSpIndex += nsp;
53  size_t nPhases = kinPtr->nPhases();
54  vector_int pLocTmp(nPhases);
55  size_t imatch = npos;
56  for (size_t ip = 0; ip < nPhases; ip++) {
57  if (ip != ns) {
58  ThermoPhase* thPtr = & kinPtr->thermo(ip);
59  if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
60  m_bulkPhases.push_back(thPtr);
61  nsp = thPtr->nSpecies();
62  m_numTotalBulkSpecies += nsp;
63  imatch = m_bulkPhases.size() - 1;
64  }
65  pLocTmp[ip] = int(imatch);
66  } else {
67  pLocTmp[ip] = -int(n);
68  }
69  }
70  pLocVec.push_back(pLocTmp);
71  }
72  m_numTotalSpecies = m_nv + m_numTotalBulkSpecies;
73  m_concSpecies.resize(m_numTotalSpecies, 0.0);
74  m_concSpeciesSave.resize(m_numTotalSpecies, 0.0);
75 
76  m_integ.reset(newIntegrator("CVODE"));
77 
78  // use backward differencing, with a full Jacobian computed
79  // numerically, and use a Newton linear iterator
80  m_integ->setMethod(BDF_Method);
81  m_integ->setProblemType(DENSE + NOJAC);
82  m_integ->setIterator(Newton_Iter);
83  m_work.resize(ntmax);
84 }
85 
86 int ImplicitSurfChem::checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* thPtr)
87 {
88  int retn = -1;
89  for (int i = 0; i < (int) m_vec.size(); i++) {
90  ThermoPhase* th = m_vec[i];
91  if (th == thPtr) {
92  return i;
93  }
94  }
95  return retn;
96 }
97 
98 void ImplicitSurfChem::getState(doublereal* c)
99 {
100  size_t loc = 0;
101  for (size_t n = 0; n < m_surf.size(); n++) {
102  m_surf[n]->getCoverages(c + loc);
103  loc += m_nsp[n];
104  }
105 }
106 
107 void ImplicitSurfChem::initialize(doublereal t0)
108 {
109  m_integ->setTolerances(m_rtol, m_atol);
110  m_integ->initialize(t0, *this);
111 }
112 
113 void ImplicitSurfChem::integrate(doublereal t0, doublereal t1)
114 {
115  m_integ->initialize(t0, *this);
116  m_integ->setMaxStepSize(t1 - t0);
117  m_integ->integrate(t1);
118  updateState(m_integ->solution());
119 }
120 
121 void ImplicitSurfChem::integrate0(doublereal t0, doublereal t1)
122 {
123  m_integ->integrate(t1);
124  updateState(m_integ->solution());
125 }
126 
127 void ImplicitSurfChem::updateState(doublereal* c)
128 {
129  size_t loc = 0;
130  for (size_t n = 0; n < m_surf.size(); n++) {
131  m_surf[n]->setCoverages(c + loc);
132  loc += m_nsp[n];
133  }
134 }
135 
136 void ImplicitSurfChem::eval(doublereal time, doublereal* y,
137  doublereal* ydot, doublereal* p)
138 {
139  updateState(y); // synchronize the surface state(s) with y
140  size_t loc = 0;
141  for (size_t n = 0; n < m_surf.size(); n++) {
142  double rs0 = 1.0/m_surf[n]->siteDensity();
143  m_vecKinPtrs[n]->getNetProductionRates(m_work.data());
144  size_t kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
145  double sum = 0.0;
146  for (size_t k = 1; k < m_nsp[n]; k++) {
147  ydot[k + loc] = m_work[kstart + k] * rs0 * m_surf[n]->size(k);
148  sum -= ydot[k];
149  }
150  ydot[loc] = sum;
151  loc += m_nsp[n];
152  }
153 }
154 
156  doublereal timeScaleOverride)
157 {
158  int ifunc;
159  // set bulkFunc. We assume that the bulk concentrations are constant.
160  int bulkFunc = BULK_ETCH;
161  // time scale - time over which to integrate equations
162  doublereal time_scale = timeScaleOverride;
163  if (!m_surfSolver) {
164  m_surfSolver.reset(new solveSP(this, bulkFunc));
165  // set ifunc, which sets the algorithm.
166  ifunc = SFLUX_INITIALIZE;
167  } else {
168  ifunc = SFLUX_RESIDUAL;
169  }
170 
171  // Possibly override the ifunc value
172  if (ifuncOverride >= 0) {
173  ifunc = ifuncOverride;
174  }
175 
176  // Get the specifications for the problem from the values
177  // in the ThermoPhase objects for all phases.
178  //
179  // 1) concentrations of all species in all phases, m_concSpecies[]
180  // 2) Temperature and pressure
183  ThermoPhase& tp = ik->thermo(0);
184  doublereal TKelvin = tp.temperature();
185  doublereal PGas = tp.pressure();
186 
187  // Make sure that there is a common temperature and pressure for all
188  // ThermoPhase objects belonging to the interfacial kinetics object, if it
189  // is required by the problem statement.
191  setCommonState_TP(TKelvin, PGas);
192  }
193 
194  doublereal reltol = 1.0E-6;
195  doublereal atol = 1.0E-20;
196 
197  // Install a filter for negative concentrations. One of the few ways solveSS
198  // can fail is if concentrations on input are below zero.
199  bool rset = false;
200  for (size_t k = 0; k < m_nv; k++) {
201  if (m_concSpecies[k] < 0.0) {
202  rset = true;
203  m_concSpecies[k] = 0.0;
204  }
205  }
206  if (rset) {
208  }
209 
210  m_surfSolver->m_ioflag = m_ioFlag;
211 
212  // Save the current solution
213  m_concSpeciesSave = m_concSpecies;
214 
215  int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
216  reltol, atol);
217  if (retn != 1) {
218  // reset the concentrations
219  m_concSpecies = m_concSpeciesSave;
220  setConcSpecies(m_concSpeciesSave.data());
221  ifunc = SFLUX_INITIALIZE;
222  retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
223  reltol, atol);
224  if (retn != 1) {
225  throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem",
226  "solveSP return an error condition!");
227  }
228  }
229 }
230 
231 void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
232 {
233  size_t kstart;
234  for (size_t ip = 0; ip < m_surf.size(); ip++) {
235  ThermoPhase* TP_ptr = m_surf[ip];
236  kstart = m_specStartIndex[ip];
237  TP_ptr->getConcentrations(vecConcSpecies + kstart);
238  }
239  kstart = m_nv;
240  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
241  ThermoPhase* TP_ptr = m_bulkPhases[ip];
242  TP_ptr->getConcentrations(vecConcSpecies + kstart);
243  kstart += TP_ptr->nSpecies();
244  }
245 }
246 
247 void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
248 {
249  size_t kstart;
250  for (size_t ip = 0; ip < m_surf.size(); ip++) {
251  ThermoPhase* TP_ptr = m_surf[ip];
252  kstart = m_specStartIndex[ip];
253  TP_ptr->setConcentrations(vecConcSpecies + kstart);
254  }
255  kstart = m_nv;
256  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
257  ThermoPhase* TP_ptr = m_bulkPhases[ip];
258  TP_ptr->setConcentrations(vecConcSpecies + kstart);
259  kstart += TP_ptr->nSpecies();
260  }
261 }
262 
263 void ImplicitSurfChem::setCommonState_TP(doublereal TKelvin, doublereal PresPa)
264 {
265  for (size_t ip = 0; ip < m_surf.size(); ip++) {
266  ThermoPhase* TP_ptr = m_surf[ip];
267  TP_ptr->setState_TP(TKelvin, PresPa);
268  }
269  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
270  ThermoPhase* TP_ptr = m_bulkPhases[ip];
271  TP_ptr->setState_TP(TKelvin, PresPa);
272  }
273 }
274 
275 }
Backward Differentiation.
Definition: Integrator.h:33
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase, assuming an ideal solution model (see Thermodynamic Properties and class SurfPhase).
std::vector< size_t > m_surfindex
index of the surface phase in each InterfaceKinetics object
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
virtual void getState(doublereal *y)
Get the current state of the solution vector.
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism...
Definition: Kinetics.h:227
std::unique_ptr< Integrator > m_integ
Pointer to the CVODE integrator.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, doublereal timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
std::vector< ThermoPhase * > m_bulkPhases
Vector of pointers to bulk phases.
Method to solve a pseudo steady state surface problem.
Definition: solveSP.h:138
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
std::unique_ptr< solveSP > m_surfSolver
Pointer to the helper method, Placid, which solves the surface problem.
std::vector< size_t > m_nsp
Vector of number of species in each Surface Phase.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition: solveSP.h:56
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:167
std::vector< SurfPhase * > m_surf
vector of pointers to surface phases.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:142
int m_ioFlag
Controls the amount of printing from this routine and underlying routines.
bool m_commonTempPressForPhases
If true, a common temperature and pressure for all surface and bulk phases associated with the surfac...
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one...
Definition: solveSP.h:21
A kinetics manager for heterogeneous reaction mechanisms.
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species...
Definition: solveSP.h:28
size_t m_nv
Total number of surface species in all surface phases.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void eval(doublereal t, doublereal *y, doublereal *ydot, doublereal *p)
Evaluate the value of ydot[k] at the current conditions.
virtual void setConcentrations(const doublereal *const conc)
Set the concentrations to the specified values within the phase.
Definition: Phase.cpp:524
void integrate0(doublereal t0, doublereal t1)
Integrate from t0 to t1 without reinitializing the integrator.
void setConcSpecies(const doublereal *const vecConcSpecies)
Sets the concentrations within phases that are unknowns in the surface problem.
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
virtual void setState_TP(doublereal t, doublereal p)
Set the temperature (K) and pressure (Pa)
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:245
void getConcSpecies(doublereal *const vecConcSpecies) const
virtual void initialize(doublereal t0=0.0)
void setCommonState_TP(doublereal TKelvin, doublereal PresPa)
Sets the state variable in all thermodynamic phases (surface and surrounding bulk phases) to the inpu...
Newton Iteration.
Definition: Integrator.h:43
void updateState(doublereal *y)
Set the mixture to a state consistent with solution vector y.
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).
std::vector< InterfaceKinetics * > m_vecKinPtrs
vector of pointers to InterfaceKinetics objects
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
void getConcentrations(doublereal *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:519
vector_fp m_concSpecies
Temporary vector - length num species in the Kinetics object.
void integrate(doublereal t0, doublereal t1)
Integrate from t0 to t1. The integrator is reinitialized first.