Cantera  3.1.0a1
Loading...
Searching...
No Matches
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
18namespace 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 */
59{
60public:
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
221protected:
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
297private:
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.
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.
vector< InterfaceKinetics * > & getObjects()
Returns a reference to the vector of pointers to the InterfaceKinetics objects.
void solvePseudoSteadyStateProblem(int ifuncOverride=-1, double timeScaleOverride=1.0)
Solve for the pseudo steady-state of the surface problem.
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.
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).