Cantera  2.5.1
ImplicitSurfChem.h
Go to the documentation of this file.
1 /**
2  * @file ImplicitSurfChem.h
3  * Declarations 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 
11 #ifndef CT_IMPSURFCHEM_H
12 #define CT_IMPSURFCHEM_H
13 
17 
18 namespace Cantera
19 {
20 
21 //! Advances the surface coverages of the associated set of SurfacePhase
22 //! objects in time
23 /*!
24  * This function advances a set of SurfPhase objects, each associated with one
25  * InterfaceKinetics object, in time. The following equation is used for each
26  * surface phase, *i*.
27  *
28  * \f[
29  * \dot \theta_k = \dot s_k (\sigma_k / s_0)
30  * \f]
31  *
32  * In this equation,
33  * - \f$ \theta_k \f$ is the site coverage for the kth species.
34  * - \f$ \dot s_k \f$ is the source term for the kth species
35  * - \f$ \sigma_k \f$ is the number of surface sites covered by each species k.
36  * - \f$ s_0 \f$ is the total site density of the interfacial phase.
37  *
38  * Additionally, the 0'th equation in the set is discarded. Instead the
39  * alternate equation is solved for
40  *
41  * \f[
42  * \sum_{k=0}^{N-1} \dot \theta_k = 0
43  * \f]
44  *
45  * This last equation serves to ensure that sum of the \f$ \theta_k \f$ values
46  * stays constant.
47  *
48  * The object uses the CVODE software to advance the surface equations.
49  *
50  * The solution vector used by this object is as follows: For each surface
51  * phase with \f$ N_s \f$ surface sites, it consists of the surface coverages
52  * \f$ \theta_k \f$ for \f$ k = 0, N_s - 1 \f$
53  *
54  * @ingroup kineticsmgr
55  */
56 class ImplicitSurfChem : public FuncEval
57 {
58 public:
59  //! Constructor for multiple surfaces.
60  /*!
61  * @param k Vector of pointers to InterfaceKinetics objects Each object
62  * consists of a surface or an edge containing internal degrees of
63  * freedom representing the concentration of surface adsorbates.
64  * @param rtol The relative tolerance for the integrator
65  * @param atol The absolute tolerance for the integrator
66  * @param maxStepSize The maximum step-size the integrator is allowed to take.
67  * If zero, this option is disabled.
68  * @param maxSteps The maximum number of time-steps the integrator can take.
69  * If not supplied, uses the default value in the CVodesIntegrator (20000).
70  * @param maxErrTestFails the maximum permissible number of error test failures
71  * If not supplied, uses the default value in CVODES (7).
72  */
73  ImplicitSurfChem(std::vector<InterfaceKinetics*> k,
74  double rtol=1.e-7, double atol=1.e-14,
75  double maxStepSize=0, size_t maxSteps=20000,
76  size_t maxErrTestFails=7);
77 
78  virtual ~ImplicitSurfChem() {};
79 
80  /*!
81  * Must be called before calling method 'advance'
82  */
83  virtual void initialize(doublereal t0 = 0.0);
84 
85  /*!
86  * Set the maximum integration step-size. Note, setting this value to zero
87  * disables this option
88  */
89  virtual void setMaxStepSize(double maxstep = 0.0);
90 
91  /*!
92  * Set the relative and absolute integration tolerances.
93  */
94  virtual void setTolerances(double rtol=1.e-7, double atol=1.e-14);
95 
96  /*!
97  * Set the maximum number of CVODES integration steps.
98  */
99  virtual void setMaxSteps(size_t maxsteps = 20000);
100 
101  /*!
102  * Set the maximum number of CVODES error test failures
103  */
104  virtual void setMaxErrTestFails(size_t maxErrTestFails = 7);
105 
106  //! Integrate from t0 to t1. The integrator is reinitialized first.
107  /*!
108  * This routine does a time accurate solve from t = t0 to t = t1.
109  * of the surface problem.
110  *
111  * @param t0 Initial Time -> this is an input
112  * @param t1 Final Time -> This is an input
113  */
114  void integrate(doublereal t0, doublereal t1);
115 
116  //! Integrate from t0 to t1 without reinitializing the integrator.
117  /*!
118  * Use when the coverages have not changed from their values on return
119  * from the last call to integrate or integrate0.
120  *
121  * @param t0 Initial Time -> this is an input
122  * @param t1 Final Time -> This is an input
123  */
124  void integrate0(doublereal t0, doublereal t1);
125 
126  //! Solve for the pseudo steady-state of the surface problem
127  /*!
128  * Solve for the steady state of the surface problem.
129  * This is the same thing as the advanceCoverages() function,
130  * but at infinite times.
131  *
132  * Note, a direct solve is carried out under the hood here,
133  * to reduce the computational time.
134  *
135  * @param ifuncOverride One of the values defined in @ref solvesp_methods.
136  * The default is -1, which means that the program will decide.
137  * @param timeScaleOverride When a pseudo transient is selected this value
138  * can be used to override the default time scale for
139  * integration which is one. When SFLUX_TRANSIENT is used, this
140  * is equal to the time over which the equations are integrated.
141  * When SFLUX_INITIALIZE is used, this is equal to the time used
142  * in the initial transient algorithm, before the equation
143  * system is solved directly.
144  */
145  void solvePseudoSteadyStateProblem(int ifuncOverride = -1,
146  doublereal timeScaleOverride = 1.0);
147 
148  // overloaded methods of class FuncEval
149 
150  //! Return the number of equations
151  virtual size_t neq() {
152  return m_nv;
153  }
154 
155  //! Evaluate the value of ydot[k] at the current conditions
156  /*!
157  * @param t Time (seconds)
158  * @param y Vector containing the current solution vector
159  * @param ydot Output vector containing the value of the
160  * derivative of the surface coverages.
161  * @param p Unused parameter pass-through parameter vector
162  */
163  virtual void eval(doublereal t, doublereal* y, doublereal* ydot,
164  doublereal* p);
165 
166  //! Get the current state of the solution vector
167  /*!
168  * @param y Value of the solution vector to be used.
169  * On output, this contains the initial value
170  * of the solution.
171  */
172  virtual void getState(doublereal* y);
173 
174  /*!
175  * Get the specifications for the problem from the values
176  * in the ThermoPhase objects for all phases.
177  *
178  * 1. concentrations of all species in all phases, #m_concSpecies
179  * 2. Temperature and pressure
180  *
181  * @param vecConcSpecies Vector of concentrations. The phase concentration
182  * vectors are contiguous within the object, in the same
183  * order as the unknown vector.
184  */
185  void getConcSpecies(doublereal* const vecConcSpecies) const;
186 
187  //! Sets the concentrations within phases that are unknowns in
188  //! the surface problem
189  /*!
190  * Fills the local concentration vector for all of the species in all of
191  * the phases that are unknowns in the surface problem.
192  *
193  * @param vecConcSpecies Vector of concentrations. The phase concentration
194  * vectors are contiguous within the object, in the same
195  * order as the unknown vector.
196  */
197  void setConcSpecies(const doublereal* const vecConcSpecies);
198 
199  //! Sets the state variable in all thermodynamic phases (surface and
200  //! surrounding bulk phases) to the input temperature and pressure
201  /*!
202  * @param TKelvin input temperature (kelvin)
203  * @param PresPa input pressure in pascal.
204  */
205  void setCommonState_TP(doublereal TKelvin, doublereal PresPa);
206 
207  //! Returns a reference to the vector of pointers to the
208  //! InterfaceKinetics objects
209  /*!
210  * This should probably go away in the future, as it opens up the class.
211  */
212  std::vector<InterfaceKinetics*> & getObjects() {
213  return m_vecKinPtrs;
214  }
215 
216  int checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* thPtr);
217 
218  void setIOFlag(int ioFlag) {
219  m_ioFlag = ioFlag;
220  }
221 
222 protected:
223  //! Set the mixture to a state consistent with solution
224  //! vector y.
225  /*!
226  * This function will set the surface site factions in the underlying
227  * SurfPhase objects to the current value of the solution vector.
228  *
229  * @param y Current value of the solution vector. The length is equal to
230  * the sum of the number of surface sites in all the surface phases.
231  */
232  void updateState(doublereal* y);
233 
234  //! vector of pointers to surface phases.
235  std::vector<SurfPhase*> m_surf;
236 
237  //! Vector of pointers to bulk phases
238  std::vector<ThermoPhase*> m_bulkPhases;
239 
240  //! vector of pointers to InterfaceKinetics objects
241  std::vector<InterfaceKinetics*> m_vecKinPtrs;
242 
243  //! Vector of number of species in each Surface Phase
244  std::vector<size_t> m_nsp;
245 
246  //! index of the surface phase in each InterfaceKinetics object
247  std::vector<size_t> m_surfindex;
248 
249  std::vector<size_t> m_specStartIndex;
250 
251  //! Total number of surface species in all surface phases
252  /*!
253  * This is the total number of unknowns in m_mode 0 problem
254  */
255  size_t m_nv;
256 
257  size_t m_numTotalBulkSpecies;
258  size_t m_numTotalSpecies;
259 
260  std::vector<vector_int> pLocVec;
261  //! Pointer to the CVODE integrator
262  std::unique_ptr<Integrator> m_integ;
263  doublereal m_atol, m_rtol; // tolerances
264  doublereal m_maxstep; //!< max step size
265  size_t m_nmax; //!< maximum number of steps allowed
266  size_t m_maxErrTestFails; //!< maximum number of error test failures allowed
267  vector_fp m_work;
268 
269  /**
270  * Temporary vector - length num species in the Kinetics object. This is
271  * the sum of the number of species in each phase included in the kinetics
272  * object.
273  */
275  vector_fp m_concSpeciesSave;
276 
277  /**
278  * Index into the species vector of the kinetics manager,
279  * pointing to the first species from the surrounding medium.
280  */
282  /**
283  * Index into the species vector of the kinetics manager, pointing to the
284  * first species from the condensed phase of the particles.
285  */
287  /**
288  * Index into the species vector of the kinetics manager, pointing to the
289  * first species from the surface of the particles
290  */
292  /**
293  * Pointer to the helper method, Placid, which solves the surface problem.
294  */
295  std::unique_ptr<solveSP> m_surfSolver;
296 
297  //! If true, a common temperature and pressure for all surface and bulk
298  //! phases associated with the surface problem is imposed
300 
301  //! We make the solveSS class a friend because we need to access all of
302  //! the above information directly. Adding the members into the class is
303  //! also a possibility.
304  friend class solveSS;
305 
306 private:
307  //! Controls the amount of printing from this routine
308  //! and underlying routines.
309  int m_ioFlag;
310 };
311 }
312 
313 #endif
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:27
Advances the surface coverages of the associated set of SurfacePhase objects in time.
void integrate0(doublereal t0, doublereal t1)
Integrate from t0 to t1 without reinitializing the integrator.
void integrate(doublereal t0, doublereal t1)
Integrate from t0 to t1. The integrator is reinitialized first.
int m_bulkSpeciesStart
Index into the species vector of the kinetics manager, pointing to the first species from the condens...
int m_mediumSpeciesStart
Index into the species vector of the kinetics manager, pointing to the first species from the surroun...
void updateState(doublereal *y)
Set the mixture to a state consistent with solution vector y.
std::vector< ThermoPhase * > m_bulkPhases
Vector of pointers to bulk phases.
void setConcSpecies(const doublereal *const vecConcSpecies)
Sets the concentrations within phases that are unknowns in the surface problem.
size_t m_nv
Total number of surface species in all surface phases.
doublereal m_maxstep
max step size
int m_ioFlag
Controls the amount of printing from this routine and underlying routines.
virtual void getState(doublereal *y)
Get the current state of the solution vector.
virtual size_t neq()
Return the number of equations.
std::vector< InterfaceKinetics * > & getObjects()
Returns a reference to the vector of pointers to the InterfaceKinetics objects.
std::unique_ptr< solveSP > m_surfSolver
Pointer to the helper method, Placid, which solves the surface problem.
std::vector< size_t > m_surfindex
index of the surface phase in each InterfaceKinetics object
bool m_commonTempPressForPhases
If true, a common temperature and pressure for all surface and bulk phases associated with the surfac...
friend class solveSS
We make the solveSS class a friend because we need to access all of the above information directly.
virtual void setTolerances(double rtol=1.e-7, double atol=1.e-14)
virtual void setMaxStepSize(double maxstep=0.0)
ImplicitSurfChem(std::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.
std::vector< SurfPhase * > m_surf
vector of pointers to surface phases.
virtual void setMaxSteps(size_t maxsteps=20000)
vector_fp m_concSpecies
Temporary vector - length num species in the Kinetics object.
void getConcSpecies(doublereal *const vecConcSpecies) const
virtual void initialize(doublereal t0=0.0)
size_t m_maxErrTestFails
maximum number of error test failures allowed
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, doublereal timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
int m_surfSpeciesStart
Index into the species vector of the kinetics manager, pointing to the first species from the surface...
std::vector< size_t > m_nsp
Vector of number of species in each Surface Phase.
virtual void setMaxErrTestFails(size_t maxErrTestFails=7)
void setCommonState_TP(doublereal TKelvin, doublereal PresPa)
Sets the state variable in all thermodynamic phases (surface and surrounding bulk phases) to the inpu...
virtual void eval(doublereal t, doublereal *y, doublereal *ydot, doublereal *p)
Evaluate the value of ydot[k] at the current conditions.
std::unique_ptr< Integrator > m_integ
Pointer to the CVODE integrator.
size_t m_nmax
maximum number of steps allowed
std::vector< InterfaceKinetics * > m_vecKinPtrs
vector of pointers to InterfaceKinetics objects
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).