Cantera  3.1.0a1
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 https://cantera.org/license.txt for license and copyright information.
10 
14 
15 namespace Cantera
16 {
17 
19  vector<InterfaceKinetics*> k, double rtol, double atol,
20  double maxStepSize, size_t maxSteps,
21  size_t maxErrTestFails) :
22  m_atol(atol),
23  m_rtol(rtol),
24  m_maxstep(maxStepSize),
25  m_nmax(maxSteps),
26  m_maxErrTestFails(maxErrTestFails)
27 {
28  size_t ntmax = 0;
29  size_t kinSpIndex = 0;
30  // Loop over the number of surface kinetics objects
31  for (size_t n = 0; n < k.size(); n++) {
32  InterfaceKinetics* kinPtr = k[n];
33  m_vecKinPtrs.push_back(kinPtr);
34  SurfPhase* surf = dynamic_cast<SurfPhase*>(&k[n]->thermo(0));
35  if (surf == nullptr) {
36  throw CanteraError("ImplicitSurfChem::ImplicitSurfChem",
37  "kinetics manager contains no surface phase");
38  }
39  m_surf.push_back(surf);
40  size_t nsp = m_surf.back()->nSpecies();
41  m_nsp.push_back(nsp);
42  m_nv += m_nsp.back();
43  size_t nt = k[n]->nTotalSpecies();
44  ntmax = std::max(nt, ntmax);
45  m_specStartIndex.push_back(kinSpIndex);
46  kinSpIndex += nsp;
47  size_t nPhases = kinPtr->nPhases();
48  vector<int> pLocTmp(nPhases);
49  pLocTmp[0] = -int(n);
50  size_t imatch = npos;
51  for (size_t ip = 1; ip < nPhases; ip++) {
52  ThermoPhase* thPtr = & kinPtr->thermo(ip);
53  if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
54  m_bulkPhases.push_back(thPtr);
55  nsp = thPtr->nSpecies();
56  m_numTotalBulkSpecies += nsp;
57  imatch = m_bulkPhases.size() - 1;
58  }
59  pLocTmp[ip] = int(imatch);
60  }
61  pLocVec.push_back(pLocTmp);
62  }
63  m_numTotalSpecies = m_nv + m_numTotalBulkSpecies;
64  m_concSpecies.resize(m_numTotalSpecies, 0.0);
65  m_concSpeciesSave.resize(m_numTotalSpecies, 0.0);
66 
67  m_integ.reset(newIntegrator("CVODE"));
68 
69  // use backward differencing, with a full Jacobian computed
70  // numerically, and use a Newton linear iterator
71  m_integ->setMethod(BDF_Method);
72  m_integ->setLinearSolverType("DENSE");
73  m_work.resize(ntmax);
74 }
75 
76 int ImplicitSurfChem::checkMatch(vector<ThermoPhase*> m_vec, ThermoPhase* thPtr)
77 {
78  int retn = -1;
79  for (int i = 0; i < (int) m_vec.size(); i++) {
80  ThermoPhase* th = m_vec[i];
81  if (th == thPtr) {
82  return i;
83  }
84  }
85  return retn;
86 }
87 
89 {
90  size_t loc = 0;
91  for (size_t n = 0; n < m_surf.size(); n++) {
92  m_surf[n]->getCoverages(c + loc);
93  loc += m_nsp[n];
94  }
95 }
96 
98 {
99  m_maxstep = maxstep;
100  if (m_maxstep > 0) {
101  m_integ->setMaxStepSize(m_maxstep);
102  }
103 }
104 
105 void ImplicitSurfChem::setTolerances(double rtol, double atol)
106 {
107  m_rtol = rtol;
108  m_atol = atol;
109  m_integ->setTolerances(m_rtol, m_atol);
110 }
111 
112 void ImplicitSurfChem::setMaxSteps(size_t maxsteps)
113 {
114  m_nmax = maxsteps;
115  m_integ->setMaxSteps(static_cast<int>(m_nmax));
116 }
117 
118 void ImplicitSurfChem::setMaxErrTestFails(size_t maxErrTestFails)
119 {
120  m_maxErrTestFails = maxErrTestFails;
121  m_integ->setMaxErrTestFails(static_cast<int>(m_maxErrTestFails));
122 }
123 
125 {
126  this->setTolerances(m_rtol, m_atol);
127  this->setMaxStepSize(m_maxstep);
128  this->setMaxSteps(m_nmax);
130  m_integ->initialize(t0, *this);
131 }
132 
133 void ImplicitSurfChem::integrate(double t0, double t1)
134 {
135  this->initialize(t0);
136  if (fabs(t1 - t0) < m_maxstep || m_maxstep == 0) {
137  // limit max step size on this run to t1 - t0
138  m_integ->setMaxStepSize(t1 - t0);
139  }
140  m_integ->integrate(t1);
141  updateState(m_integ->solution());
142 }
143 
144 void ImplicitSurfChem::integrate0(double t0, double t1)
145 {
146  m_integ->integrate(t1);
147  updateState(m_integ->solution());
148 }
149 
151 {
152  size_t loc = 0;
153  for (size_t n = 0; n < m_surf.size(); n++) {
154  m_surf[n]->setCoverages(c + loc);
155  loc += m_nsp[n];
156  }
157 }
158 
159 void ImplicitSurfChem::eval(double time, double* y, double* ydot, double* p)
160 {
161  updateState(y); // synchronize the surface state(s) with y
162  size_t loc = 0;
163  for (size_t n = 0; n < m_surf.size(); n++) {
164  double rs0 = 1.0/m_surf[n]->siteDensity();
165  m_vecKinPtrs[n]->getNetProductionRates(m_work.data());
166  double sum = 0.0;
167  for (size_t k = 1; k < m_nsp[n]; k++) {
168  ydot[k + loc] = m_work[k] * rs0 * m_surf[n]->size(k);
169  sum -= ydot[k];
170  }
171  ydot[loc] = sum;
172  loc += m_nsp[n];
173  }
174 }
175 
177  double timeScaleOverride)
178 {
179  int ifunc;
180  // set bulkFunc. We assume that the bulk concentrations are constant.
181  int bulkFunc = BULK_ETCH;
182  // time scale - time over which to integrate equations
183  double time_scale = timeScaleOverride;
184  if (!m_surfSolver) {
185  m_surfSolver = make_unique<solveSP>(this, bulkFunc);
186  // set ifunc, which sets the algorithm.
187  ifunc = SFLUX_INITIALIZE;
188  } else {
189  ifunc = SFLUX_RESIDUAL;
190  }
191 
192  // Possibly override the ifunc value
193  if (ifuncOverride >= 0) {
194  ifunc = ifuncOverride;
195  }
196 
197  // Get the specifications for the problem from the values
198  // in the ThermoPhase objects for all phases.
199  //
200  // 1) concentrations of all species in all phases, m_concSpecies[]
201  // 2) Temperature and pressure
204  ThermoPhase& tp = ik->thermo(0);
205  double TKelvin = tp.temperature();
206  double PGas = tp.pressure();
207 
208  // Make sure that there is a common temperature and pressure for all
209  // ThermoPhase objects belonging to the interfacial kinetics object, if it
210  // is required by the problem statement.
212  setCommonState_TP(TKelvin, PGas);
213  }
214 
215  double reltol = 1.0E-6;
216  double atol = 1.0E-20;
217 
218  // Install a filter for negative concentrations. One of the few ways solveSP
219  // can fail is if concentrations on input are below zero.
220  bool rset = false;
221  for (size_t k = 0; k < m_nv; k++) {
222  if (m_concSpecies[k] < 0.0) {
223  rset = true;
224  m_concSpecies[k] = 0.0;
225  }
226  }
227  if (rset) {
229  }
230 
231  m_surfSolver->m_ioflag = m_ioFlag;
232 
233  // Save the current solution
234  m_concSpeciesSave = m_concSpecies;
235 
236  int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
237  reltol, atol);
238  if (retn != 1) {
239  // reset the concentrations
240  m_concSpecies = m_concSpeciesSave;
241  setConcSpecies(m_concSpeciesSave.data());
242  ifunc = SFLUX_INITIALIZE;
243  retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
244  reltol, atol);
245  if (retn != 1) {
246  throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem",
247  "solveSP return an error condition!");
248  }
249  }
250 }
251 
252 void ImplicitSurfChem::getConcSpecies(double* const vecConcSpecies) const
253 {
254  size_t kstart;
255  for (size_t ip = 0; ip < m_surf.size(); ip++) {
256  ThermoPhase* TP_ptr = m_surf[ip];
257  kstart = m_specStartIndex[ip];
258  TP_ptr->getConcentrations(vecConcSpecies + kstart);
259  }
260  kstart = m_nv;
261  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
262  ThermoPhase* TP_ptr = m_bulkPhases[ip];
263  TP_ptr->getConcentrations(vecConcSpecies + kstart);
264  kstart += TP_ptr->nSpecies();
265  }
266 }
267 
268 void ImplicitSurfChem::setConcSpecies(const double* const vecConcSpecies)
269 {
270  size_t kstart;
271  for (size_t ip = 0; ip < m_surf.size(); ip++) {
272  ThermoPhase* TP_ptr = m_surf[ip];
273  kstart = m_specStartIndex[ip];
274  TP_ptr->setConcentrations(vecConcSpecies + kstart);
275  }
276  kstart = m_nv;
277  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
278  ThermoPhase* TP_ptr = m_bulkPhases[ip];
279  TP_ptr->setConcentrations(vecConcSpecies + kstart);
280  kstart += TP_ptr->nSpecies();
281  }
282 }
283 
284 void ImplicitSurfChem::setCommonState_TP(double TKelvin, double PresPa)
285 {
286  for (size_t ip = 0; ip < m_surf.size(); ip++) {
287  ThermoPhase* TP_ptr = m_surf[ip];
288  TP_ptr->setState_TP(TKelvin, PresPa);
289  }
290  for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
291  ThermoPhase* TP_ptr = m_bulkPhases[ip];
292  TP_ptr->setState_TP(TKelvin, PresPa);
293  }
294 }
295 
296 }
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,...
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
vector< ThermoPhase * > m_bulkPhases
Vector of pointers to bulk phases.
void eval(double t, double *y, double *ydot, double *p) override
Evaluate the value of ydot[k] at the current conditions.
void setConcSpecies(const double *const vecConcSpecies)
Sets the concentrations within phases that are unknowns in the surface problem.
void setCommonState_TP(double TKelvin, double PresPa)
Sets the state variable in all thermodynamic phases (surface and surrounding bulk phases) to the inpu...
size_t m_nv
Total number of surface species in all surface phases.
int m_ioFlag
Controls the amount of printing from this routine and underlying routines.
void initialize(double t0=0.0)
Must be called before calling method 'advance'.
void integrate(double t0, double t1)
Integrate from t0 to t1. The integrator is reinitialized first.
void getConcSpecies(double *const vecConcSpecies) const
Get the specifications for the problem from the values in the ThermoPhase objects for all phases.
vector< SurfPhase * > m_surf
vector of pointers to surface phases.
vector< InterfaceKinetics * > m_vecKinPtrs
vector of pointers to InterfaceKinetics objects
unique_ptr< solveSP > m_surfSolver
Pointer to the helper method, Placid, which solves the surface problem.
double m_maxstep
max step size
bool m_commonTempPressForPhases
If true, a common temperature and pressure for all surface and bulk phases associated with the surfac...
void setTolerances(double rtol=1.e-7, double atol=1.e-14)
Set the relative and absolute integration tolerances.
vector< size_t > m_nsp
Vector of number of species in each Surface Phase.
void setMaxStepSize(double maxstep=0.0)
Set the maximum integration step-size.
vector< double > m_concSpecies
Temporary vector - length num species in the Kinetics object.
void setMaxSteps(size_t maxsteps=20000)
Set the maximum number of CVODES integration steps.
void updateState(double *y)
Set the mixture to a state consistent with solution vector y.
void getState(double *y) override
Get the current state of the solution vector.
size_t m_maxErrTestFails
maximum number of error test failures allowed
unique_ptr< Integrator > m_integ
Pointer to the CVODE integrator.
void integrate0(double t0, double t1)
Integrate from t0 to t1 without reinitializing the integrator.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, double timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
ImplicitSurfChem(vector< InterfaceKinetics * > k, double rtol=1.e-7, double atol=1.e-14, double maxStepSize=0, size_t maxSteps=20000, size_t maxErrTestFails=7)
Constructor for multiple surfaces.
void setMaxErrTestFails(size_t maxErrTestFails=7)
Set the maximum number of CVODES error test failures.
size_t m_nmax
maximum number of steps allowed
A kinetics manager for heterogeneous reaction mechanisms.
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:184
ThermoPhase & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition: Kinetics.h:242
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:482
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
double temperature() const
Temperature (K).
Definition: Phase.h:562
virtual void setConcentrations(const double *const conc)
Set the concentrations to the specified values within the phase.
Definition: Phase.cpp:487
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition: Phase.h:580
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:98
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
virtual void setState_TP(double t, double p)
Set the temperature (K) and pressure (Pa)
Integrator * newIntegrator(const string &itype)
Create new Integrator object.
Definition: Integrators.cpp:14
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition: solveSP.h:58
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_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
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
@ BDF_Method
Backward Differentiation.
Definition: Integrator.h:24
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).