Cantera  3.1.0a1
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 //! @defgroup surfSolverGroup Surface Problem Solver
22 
23 //! Advances the surface coverages of the associated set of SurfacePhase
24 //! objects in time
25 /*!
26  * This function advances a set of SurfPhase objects, each associated with one
27  * InterfaceKinetics object, in time. The following equation is used for each
28  * surface phase, *i*.
29  *
30  * @f[
31  * \dot \theta_k = \dot s_k (\sigma_k / s_0)
32  * @f]
33  *
34  * In this equation,
35  * - @f$ \theta_k @f$ is the site coverage for the kth species.
36  * - @f$ \dot s_k @f$ is the source term for the kth species
37  * - @f$ \sigma_k @f$ is the number of surface sites covered by each species k.
38  * - @f$ s_0 @f$ is the total site density of the interfacial phase.
39  *
40  * Additionally, the 0'th equation in the set is discarded. Instead the
41  * alternate equation is solved for
42  *
43  * @f[
44  * \sum_{k=0}^{N-1} \dot \theta_k = 0
45  * @f]
46  *
47  * This last equation serves to ensure that sum of the @f$ \theta_k @f$ values
48  * stays constant.
49  *
50  * The object uses the CVODE software to advance the surface equations.
51  *
52  * The solution vector used by this object is as follows: For each surface
53  * phase with @f$ N_s @f$ surface sites, it consists of the surface coverages
54  * @f$ \theta_k @f$ for @f$ k = 0, N_s - 1 @f$
55  *
56  * @ingroup surfSolverGroup
57  */
58 class ImplicitSurfChem : public FuncEval
59 {
60 public:
61  //! Constructor for multiple surfaces.
62  /*!
63  * @param k Vector of pointers to InterfaceKinetics objects Each object
64  * consists of a surface or an edge containing internal degrees of
65  * freedom representing the concentration of surface adsorbates.
66  * @param rtol The relative tolerance for the integrator
67  * @param atol The absolute tolerance for the integrator
68  * @param maxStepSize The maximum step-size the integrator is allowed to take.
69  * If zero, this option is disabled.
70  * @param maxSteps The maximum number of time-steps the integrator can take.
71  * If not supplied, uses the default value in the CVodesIntegrator (20000).
72  * @param maxErrTestFails the maximum permissible number of error test failures
73  * If not supplied, uses the default value in CVODES (7).
74  */
75  ImplicitSurfChem(vector<InterfaceKinetics*> k,
76  double rtol=1.e-7, double atol=1.e-14,
77  double maxStepSize=0, size_t maxSteps=20000,
78  size_t maxErrTestFails=7);
79 
80  /**
81  * Must be called before calling method 'advance'
82  */
83  void initialize(double t0=0.0);
84 
85  /**
86  * Set the maximum integration step-size. Note, setting this value to zero
87  * disables this option
88  */
89  void setMaxStepSize(double maxstep=0.0);
90 
91  /**
92  * Set the relative and absolute integration tolerances.
93  */
94  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  void setMaxSteps(size_t maxsteps=20000);
100 
101  /**
102  * Set the maximum number of CVODES error test failures
103  */
104  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(double t0, double 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(double t0, double 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  double timeScaleOverride = 1.0);
147 
148  // overloaded methods of class FuncEval
149 
150  //! Return the number of equations
151  size_t neq() const override {
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  void eval(double t, double* y, double* ydot, double* p) override;
164 
165  //! Get the current state of the solution vector
166  /*!
167  * @param y Value of the solution vector to be used.
168  * On output, this contains the initial value
169  * of the solution.
170  */
171  void getState(double* y) override;
172 
173  /**
174  * Get the specifications for the problem from the values
175  * in the ThermoPhase objects for all phases.
176  *
177  * 1. concentrations of all species in all phases, #m_concSpecies
178  * 2. Temperature and pressure
179  *
180  * @param vecConcSpecies Vector of concentrations. The phase concentration
181  * vectors are contiguous within the object, in the same
182  * order as the unknown vector.
183  */
184  void getConcSpecies(double* const vecConcSpecies) const;
185 
186  //! Sets the concentrations within phases that are unknowns in
187  //! the surface problem
188  /*!
189  * Fills the local concentration vector for all of the species in all of
190  * the phases that are unknowns in the surface problem.
191  *
192  * @param vecConcSpecies Vector of concentrations. The phase concentration
193  * vectors are contiguous within the object, in the same
194  * order as the unknown vector.
195  */
196  void setConcSpecies(const double* const vecConcSpecies);
197 
198  //! Sets the state variable in all thermodynamic phases (surface and
199  //! surrounding bulk phases) to the input temperature and pressure
200  /*!
201  * @param TKelvin input temperature (kelvin)
202  * @param PresPa input pressure in pascal.
203  */
204  void setCommonState_TP(double TKelvin, double PresPa);
205 
206  //! Returns a reference to the vector of pointers to the
207  //! InterfaceKinetics objects
208  /*!
209  * This should probably go away in the future, as it opens up the class.
210  */
211  vector<InterfaceKinetics*> & getObjects() {
212  return m_vecKinPtrs;
213  }
214 
215  int checkMatch(vector<ThermoPhase*> m_vec, ThermoPhase* thPtr);
216 
217  void setIOFlag(int ioFlag) {
218  m_ioFlag = ioFlag;
219  }
220 
221 protected:
222  //! Set the mixture to a state consistent with solution
223  //! vector y.
224  /*!
225  * This function will set the surface site factions in the underlying
226  * SurfPhase objects to the current value of the solution vector.
227  *
228  * @param y Current value of the solution vector. The length is equal to
229  * the sum of the number of surface sites in all the surface phases.
230  */
231  void updateState(double* y);
232 
233  //! vector of pointers to surface phases.
234  vector<SurfPhase*> m_surf;
235 
236  //! Vector of pointers to bulk phases
237  vector<ThermoPhase*> m_bulkPhases;
238 
239  //! vector of pointers to InterfaceKinetics objects
240  vector<InterfaceKinetics*> m_vecKinPtrs;
241 
242  //! Vector of number of species in each Surface Phase
243  vector<size_t> m_nsp;
244 
245  vector<size_t> m_specStartIndex;
246 
247  //! Total number of surface species in all surface phases
248  /*!
249  * This is the total number of unknowns in m_mode 0 problem
250  */
251  size_t m_nv = 0;
252 
253  size_t m_numTotalBulkSpecies = 0;
254  size_t m_numTotalSpecies = 0;
255 
256  vector<vector<int>> pLocVec;
257  //! Pointer to the CVODE integrator
258  unique_ptr<Integrator> m_integ;
259  double m_atol, m_rtol; // tolerances
260  double m_maxstep; //!< max step size
261  size_t m_nmax; //!< maximum number of steps allowed
262  size_t m_maxErrTestFails; //!< maximum number of error test failures allowed
263  vector<double> m_work;
264 
265  /**
266  * Temporary vector - length num species in the Kinetics object. This is
267  * the sum of the number of species in each phase included in the kinetics
268  * object.
269  */
270  vector<double> m_concSpecies;
271  vector<double> m_concSpeciesSave;
272 
273  /**
274  * Index into the species vector of the kinetics manager,
275  * pointing to the first species from the surrounding medium.
276  */
278  /**
279  * Index into the species vector of the kinetics manager, pointing to the
280  * first species from the condensed phase of the particles.
281  */
283  /**
284  * Index into the species vector of the kinetics manager, pointing to the
285  * first species from the surface of the particles
286  */
288  /**
289  * Pointer to the helper method, Placid, which solves the surface problem.
290  */
291  unique_ptr<solveSP> m_surfSolver;
292 
293  //! If true, a common temperature and pressure for all surface and bulk
294  //! phases associated with the surface problem is imposed
296 
297 private:
298  //! Controls the amount of printing from this routine
299  //! and underlying routines.
300  int m_ioFlag = 0;
301 };
302 }
303 
304 #endif
Virtual base class for ODE/DAE right-hand-side function evaluators.
Definition: FuncEval.h:32
Advances the surface coverages of the associated set of SurfacePhase objects in time.
vector< ThermoPhase * > m_bulkPhases
Vector of pointers to bulk phases.
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 eval(double t, double *y, double *ydot, double *p) override
Evaluate the value of ydot[k] at the current conditions.
vector< InterfaceKinetics * > & getObjects()
Returns a reference to the vector of pointers to the InterfaceKinetics objects.
size_t neq() const override
Return the number of equations.
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.
int m_surfSpeciesStart
Index into the species vector of the kinetics manager, pointing to the first species from the surface...
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
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).