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 14 #include "cantera/numerics/Integrator.h" 15 #include "cantera/kinetics/InterfaceKinetics.h" 16 #include "cantera/kinetics/solveSP.h" 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 ~ImplicitSurfChem()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 neq()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 */ getObjects()212 std::vector<InterfaceKinetics*> & getObjects() { 213 return m_vecKinPtrs; 214 } 215 216 int checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* thPtr); 217 setIOFlag(int ioFlag)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 */ 274 vector_fp m_concSpecies; 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 */ 281 int m_mediumSpeciesStart; 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 */ 286 int m_bulkSpeciesStart; 287 /** 288 * Index into the species vector of the kinetics manager, pointing to the 289 * first species from the surface of the particles 290 */ 291 int m_surfSpeciesStart; 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 299 bool m_commonTempPressForPhases; 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 314