Cantera  2.1.2
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 // Copyright 2001 California Institute of Technology
8 
12 
13 using namespace std;
14 
15 namespace Cantera
16 {
17 
18 ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
19  FuncEval(),
20  m_nsurf(0),
21  m_nv(0),
22  m_numBulkPhases(0),
23  m_numTotalBulkSpecies(0),
24  m_numTotalSpecies(0),
25  m_integ(0),
26  m_atol(1.e-14),
27  m_rtol(1.e-7),
28  m_maxstep(0.0),
29  m_mediumSpeciesStart(-1),
30  m_bulkSpeciesStart(-1),
31  m_surfSpeciesStart(-1),
32  m_surfSolver(0),
33  m_commonTempPressForPhases(true),
34  m_ioFlag(0)
35 {
36  m_nsurf = k.size();
37  size_t ns, nsp;
38  size_t nt, ntmax = 0;
39  size_t kinSpIndex = 0;
40  // Loop over the number of surface kinetics objects
41  for (size_t n = 0; n < m_nsurf; n++) {
42  InterfaceKinetics* kinPtr = k[n];
43  m_vecKinPtrs.push_back(kinPtr);
44  ns = k[n]->surfacePhaseIndex();
45  if (ns == npos)
46  throw CanteraError("ImplicitSurfChem",
47  "kinetics manager contains no surface phase");
48  m_surfindex.push_back(ns);
49  m_surf.push_back((SurfPhase*)&k[n]->thermo(ns));
50  nsp = m_surf.back()->nSpecies();
51  m_nsp.push_back(nsp);
52  m_nv += m_nsp.back();
53  nt = k[n]->nTotalSpecies();
54  if (nt > ntmax) {
55  ntmax = nt;
56  }
57  m_specStartIndex.push_back(kinSpIndex);
58  kinSpIndex += nsp;
59 
60  size_t nPhases = kinPtr->nPhases();
61  vector_int pLocTmp(nPhases);
62  size_t imatch = npos;
63  for (size_t ip = 0; ip < nPhases; ip++) {
64  if (ip != ns) {
65  ThermoPhase* thPtr = & kinPtr->thermo(ip);
66  if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
67  m_bulkPhases.push_back(thPtr);
68  m_numBulkPhases++;
69  nsp = thPtr->nSpecies();
70  m_nspBulkPhases.push_back(nsp);
71  m_numTotalBulkSpecies += nsp;
72  imatch = m_bulkPhases.size() - 1;
73  }
74  pLocTmp[ip] = int(imatch);
75  } else {
76  pLocTmp[ip] = -int(n);
77  }
78  }
79  pLocVec.push_back(pLocTmp);
80 
81  }
82  m_numTotalSpecies = m_nv + m_numTotalBulkSpecies;
83  m_concSpecies.resize(m_numTotalSpecies, 0.0);
84  m_concSpeciesSave.resize(m_numTotalSpecies, 0.0);
85 
86  m_integ = newIntegrator("CVODE");
87 
88 
89 
90  // use backward differencing, with a full Jacobian computed
91  // numerically, and use a Newton linear iterator
92 
94  m_integ->setProblemType(DENSE + NOJAC);
96  m_work.resize(ntmax);
97 }
98 
99 int ImplicitSurfChem::checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* thPtr)
100 {
101  int retn = -1;
102  for (int i = 0; i < (int) m_vec.size(); i++) {
103  ThermoPhase* th = m_vec[i];
104  if (th == thPtr) {
105  return i;
106  }
107  }
108  return retn;
109 }
110 
112 {
113  delete m_integ;
114  delete m_surfSolver;
115 }
116 
117 void ImplicitSurfChem::getInitialConditions(doublereal t0, size_t lenc,
118  doublereal* c)
119 {
120  size_t loc = 0;
121  for (size_t n = 0; n < m_nsurf; n++) {
122  m_surf[n]->getCoverages(c + loc);
123  loc += m_nsp[n];
124  }
125 }
126 
127 void ImplicitSurfChem::initialize(doublereal t0)
128 {
129  m_integ->setTolerances(m_rtol, m_atol);
130  m_integ->initialize(t0, *this);
131 }
132 
133 void ImplicitSurfChem::integrate(doublereal t0, doublereal t1)
134 {
135  m_integ->initialize(t0, *this);
136  m_integ->setMaxStepSize(t1 - t0);
137  m_integ->integrate(t1);
139 }
140 
141 void ImplicitSurfChem::integrate0(doublereal t0, doublereal t1)
142 {
143  m_integ->integrate(t1);
145 }
146 
147 void ImplicitSurfChem::updateState(doublereal* c)
148 {
149  size_t loc = 0;
150  for (size_t n = 0; n < m_nsurf; n++) {
151  m_surf[n]->setCoverages(c + loc);
152  loc += m_nsp[n];
153  }
154 }
155 
156 void ImplicitSurfChem::eval(doublereal time, doublereal* y,
157  doublereal* ydot, doublereal* p)
158 {
159  updateState(y); // synchronize the surface state(s) with y
160  doublereal rs0, sum;
161  size_t loc, kstart;
162  for (size_t n = 0; n < m_nsurf; n++) {
163  rs0 = 1.0/m_surf[n]->siteDensity();
164  m_vecKinPtrs[n]->getNetProductionRates(DATA_PTR(m_work));
165  kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
166  sum = 0.0;
167  loc = 0;
168  for (size_t k = 1; k < m_nsp[n]; k++) {
169  ydot[k + loc] = m_work[kstart + k] * rs0 * m_surf[n]->size(k);
170  sum -= ydot[k];
171  }
172  ydot[loc] = sum;
173  loc += m_nsp[n];
174  }
175 }
176 
178  doublereal timeScaleOverride)
179 {
180  int ifunc;
181  /*
182  * set bulkFunc
183  * -> We assume that the bulk concentrations are constant.
184  */
185  int bulkFunc = BULK_ETCH;
186  /*
187  * time scale - time over which to integrate equations
188  */
189  doublereal time_scale = timeScaleOverride;
190  if (!m_surfSolver) {
191  m_surfSolver = new solveSP(this, bulkFunc);
192  /*
193  * set ifunc, which sets the algorithm.
194  */
195  ifunc = SFLUX_INITIALIZE;
196  } else {
197  ifunc = SFLUX_RESIDUAL;
198  }
199 
200  // Possibly override the ifunc value
201  if (ifuncOverride >= 0) {
202  ifunc = ifuncOverride;
203  }
204 
205  /*
206  * Get the specifications for the problem from the values
207  * in the ThermoPhase objects for all phases.
208  *
209  * 1) concentrations of all species in all phases, m_concSpecies[]
210  * 2) Temperature and pressure
211  */
214  ThermoPhase& tp = ik->thermo(0);
215  doublereal TKelvin = tp.temperature();
216  doublereal PGas = tp.pressure();
217  /*
218  * Make sure that there is a common temperature and
219  * pressure for all ThermoPhase objects belonging to the
220  * interfacial kinetics object, if it is required by
221  * the problem statement.
222  */
224  setCommonState_TP(TKelvin, PGas);
225  }
226 
227  doublereal reltol = 1.0E-6;
228  doublereal atol = 1.0E-20;
229 
230  /*
231  * Install a filter for negative concentrations. One of the
232  * few ways solveSS can fail is if concentrations on input
233  * are below zero.
234  */
235  bool rset = false;
236  for (size_t k = 0; k < m_nv; k++) {
237  if (m_concSpecies[k] < 0.0) {
238  rset = true;
239  m_concSpecies[k] = 0.0;
240  }
241  }
242  if (rset) {
244  }
245 
246  m_surfSolver->m_ioflag = m_ioFlag;
247 
248  // Save the current solution
249  copy(m_concSpecies.begin(), m_concSpecies.end(), m_concSpeciesSave.begin());
250 
251 
252  int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
253  reltol, atol);
254  if (retn != 1) {
255  // reset the concentrations
256  copy(m_concSpeciesSave.begin(), m_concSpeciesSave.end(), m_concSpecies.begin());
257  setConcSpecies(DATA_PTR(m_concSpeciesSave));
258  ifunc = SFLUX_INITIALIZE;
259  retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
260  reltol, atol);
261 
262  if (retn != 1) {
263  throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem",
264  "solveSP return an error condition!");
265  }
266  }
267 }
268 
269 void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
270 {
271  size_t kstart;
272  for (size_t ip = 0; ip < m_nsurf; ip++) {
273  ThermoPhase* TP_ptr = m_surf[ip];
274  kstart = m_specStartIndex[ip];
275  TP_ptr->getConcentrations(vecConcSpecies + kstart);
276  }
277  kstart = m_nv;
278  for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
279  ThermoPhase* TP_ptr = m_bulkPhases[ip];
280  TP_ptr->getConcentrations(vecConcSpecies + kstart);
281  kstart += TP_ptr->nSpecies();
282  }
283 }
284 
285 void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
286 {
287  size_t kstart;
288  for (size_t ip = 0; ip < m_nsurf; ip++) {
289  ThermoPhase* TP_ptr = m_surf[ip];
290  kstart = m_specStartIndex[ip];
291  TP_ptr->setConcentrations(vecConcSpecies + kstart);
292  }
293  kstart = m_nv;
294  for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
295  ThermoPhase* TP_ptr = m_bulkPhases[ip];
296  TP_ptr->setConcentrations(vecConcSpecies + kstart);
297  kstart += TP_ptr->nSpecies();
298  }
299 }
300 
302 setCommonState_TP(doublereal TKelvin, doublereal PresPa)
303 {
304  for (size_t ip = 0; ip < m_nsurf; ip++) {
305  ThermoPhase* TP_ptr = m_surf[ip];
306  TP_ptr->setState_TP(TKelvin, PresPa);
307  }
308  for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
309  ThermoPhase* TP_ptr = m_bulkPhases[ip];
310  TP_ptr->setState_TP(TKelvin, PresPa);
311  }
312 }
313 
314 }
Backward Differentiation.
Definition: Integrator.h:33
Integrator * m_integ
Pointer to the cvode integrator.
std::vector< size_t > m_surfindex
index of the surface phase in each InterfaceKinetics object
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:288
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, doublereal timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
Definition: Integrator.h:176
virtual void integrate(doublereal tout)
Integrate the system of equations.
Definition: Integrator.h:122
std::vector< ThermoPhase * > m_bulkPhases
Vector of pointers to bulk phases.
Method to solve a pseudo steady state surface problem.
Definition: solveSP.h:145
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
std::vector< size_t > m_nsp
Vector of number of species in each Surface Phase.
virtual void setTolerances(doublereal reltol, size_t n, doublereal *abstol)
Set or reset the number of equations.
Definition: Integrator.h:73
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition: solveSP.h:60
void getConcentrations(doublereal *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:572
virtual void initialize(doublereal t0, FuncEval &func)
Initialize the integrator for a new problem.
Definition: Integrator.h:109
solveSP * m_surfSolver
Pointer to the helper method, Placid, which solves the surface problem.
std::vector< SurfPhase * > m_surf
vector of pointers to surface phases.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:167
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:143
int solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, doublereal PGas, doublereal reltol, doublereal abstol)
Main routine that actually calculates the pseudo steady state of the surface problem.
Definition: solveSP.cpp:141
int m_ioFlag
Controls the amount of printing from this routine and underlying routines.
virtual void setProblemType(int probtype)
Set the problem type.
Definition: Integrator.h:99
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:25
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:228
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:32
virtual ~ImplicitSurfChem()
Destructor.
virtual void getInitialConditions(doublereal t0, size_t leny, doublereal *y)
Set the initial conditions for the solution vector.
size_t m_nv
Total number of surface species in all surface phases.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
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:577
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)
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:314
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
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...
virtual doublereal & solution(size_t k)
The current value of the solution of equation k.
Definition: Integrator.h:137
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).
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
size_t m_nsurf
Total number of surface phases.
std::vector< InterfaceKinetics * > m_vecKinPtrs
vector of pointers to InterfaceKinetics objects
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:23
virtual void setIterator(IterType t)
Set the linear iterator.
Definition: Integrator.h:171
virtual void setMethod(MethodType t)
Set the solution method.
Definition: Integrator.h:166
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.