diff --git a/Cantera/src/kinetics/ImplicitSurfChem.cpp b/Cantera/src/kinetics/ImplicitSurfChem.cpp index 2c2122b0c..98003b194 100755 --- a/Cantera/src/kinetics/ImplicitSurfChem.cpp +++ b/Cantera/src/kinetics/ImplicitSurfChem.cpp @@ -21,72 +21,135 @@ #include "ImplicitSurfChem.h" #include "Integrator.h" +#include "solveSP.h" using namespace std; namespace Cantera { + + // Constructor + ImplicitSurfChem::ImplicitSurfChem(vector k) : + FuncEval(), + m_nsurf(0), + m_nv(0), + m_numBulkPhases(0), + m_numTotalBulkSpecies(0), + m_numTotalSpecies(0), + m_integ(0), + m_atol(1.e-14), + m_rtol(1.e-7), + m_maxstep(0.0), + m_mediumSpeciesStart(-1), + m_bulkSpeciesStart(-1), + m_surfSpeciesStart(-1), + m_surfSolver(0), + m_commonTempPressForPhases(true), + m_ioFlag(0) + { + m_nsurf = static_cast(k.size()); + int ns, nsp; + int nt, ntmax = 0; + int kinSpIndex = 0; + // Loop over the number of surface kinetics objects + for (int n = 0; n < m_nsurf; n++) { + InterfaceKinetics *kinPtr = k[n]; + m_vecKinPtrs.push_back(kinPtr); + ns = k[n]->surfacePhaseIndex(); + if (ns < 0) + throw CanteraError("ImplicitSurfChem", + "kinetics manager contains no surface phase"); + m_surfindex.push_back(ns); + m_surf.push_back((SurfPhase*)&k[n]->thermo(ns)); + nsp = m_surf.back()->nSpecies(); + m_nsp.push_back(nsp); + m_nv += m_nsp.back(); + nt = k[n]->nTotalSpecies(); + if (nt > ntmax) ntmax = nt; + m_specStartIndex.push_back(kinSpIndex); + kinSpIndex += nsp; - - - - ImplicitSurfChem::ImplicitSurfChem(vector k) - : FuncEval(), m_nv(0), m_integ(0), - m_atol(1.e-14), m_rtol(1.e-7), m_maxstep(0.0) - { - m_nsurf = static_cast(k.size()); - int ns; - int nt, ntmax = 0; - for (int n = 0; n < m_nsurf; n++) { - m_kin.push_back(k[n]); - ns = k[n]->surfacePhaseIndex(); - if (ns < 0) - throw CanteraError("ImplicitSurfChem", - "kinetics manager contains no surface phase"); - m_surfindex.push_back(ns); - m_surf.push_back((SurfPhase*)&k[n]->thermo(ns)); - m_nsp.push_back(m_surf.back()->nSpecies()); - m_nv += m_nsp.back(); - nt = k[n]->nTotalSpecies(); - if (nt > ntmax) ntmax = nt; - } - m_integ = newIntegrator("CVODE");// CVodeInt; - - // use backward differencing, with a full Jacobian computed - // numerically, and use a Newton linear iterator - - m_integ->setMethod(BDF_Method); - m_integ->setProblemType(DENSE + NOJAC); - m_integ->setIterator(Newton_Iter); - m_work.resize(ntmax); + int nPhases = kinPtr->nPhases(); + vector_int pLocTmp(nPhases); + int imatch = -1; + for (int ip = 0; ip < nPhases; ip++) { + if (ip != ns) { + ThermoPhase *thPtr = & kinPtr->thermo(ip); + if ((imatch = checkMatch(m_bulkPhases, thPtr)) < 0) { + m_bulkPhases.push_back(thPtr); + m_numBulkPhases++; + nsp = thPtr->nSpecies(); + m_nspBulkPhases.push_back(nsp); + m_numTotalBulkSpecies += nsp; + imatch = m_bulkPhases.size() - 1; + } + pLocTmp[ip] = imatch; + } else { + pLocTmp[ip] = -n; + } + } + pLocVec.push_back(pLocTmp); + } + m_numTotalSpecies = m_nv + m_numTotalBulkSpecies; + m_concSpecies.resize(m_numTotalSpecies, 0.0); + m_concSpeciesSave.resize(m_numTotalSpecies, 0.0); + + m_integ = newIntegrator("CVODE"); - /** - * Destructor. Deletes the integrator. - */ - ImplicitSurfChem::~ImplicitSurfChem(){ - delete m_integ; + + + // use backward differencing, with a full Jacobian computed + // numerically, and use a Newton linear iterator + + m_integ->setMethod(BDF_Method); + m_integ->setProblemType(DENSE + NOJAC); + m_integ->setIterator(Newton_Iter); + m_work.resize(ntmax); } - // overloaded method of FuncEval. Called by the integrator to - // get the initial conditions. - void ImplicitSurfChem::getInitialConditions(double t0, size_t lenc, - double* c) - { - int loc = 0; - for (int n = 0; n < m_nsurf; n++) { - m_surf[n]->getCoverages(c + loc); - loc += m_nsp[n]; - } + int ImplicitSurfChem::checkMatch(std::vector m_vec, ThermoPhase *thPtr) { + int retn = -1; + for (int i = 0; i < (int) m_vec.size(); i++) { + ThermoPhase *th = m_vec[i]; + if (th == thPtr) { + return i; + } } + return retn; + } + + /* + * Destructor. Deletes the integrator. + */ + ImplicitSurfChem::~ImplicitSurfChem(){ + if (m_integ) { + delete m_integ; + } + if (m_surfSolver) { + delete m_surfSolver; + } + } + + // overloaded method of FuncEval. Called by the integrator to + // get the initial conditions. + void ImplicitSurfChem::getInitialConditions(double t0, size_t lenc, + double* c) + { + int loc = 0; + for (int n = 0; n < m_nsurf; n++) { + m_surf[n]->getCoverages(c + loc); + loc += m_nsp[n]; + } + } - /** - * Must be called before calling method 'advance' - */ - void ImplicitSurfChem::initialize(doublereal t0) { - m_integ->setTolerances(m_rtol, m_atol); - m_integ->initialize(t0, *this); - } + /* + * Must be called before calling method 'advance' + */ + void ImplicitSurfChem::initialize(doublereal t0) { + m_integ->setTolerances(m_rtol, m_atol); + m_integ->initialize(t0, *this); + } // Integrate from t0 to t1. The integrator is reinitialized first. /* @@ -117,16 +180,15 @@ namespace Cantera { updateState(m_integ->solution()); } - void ImplicitSurfChem::updateState(doublereal* c) { - int loc = 0; - for (int n = 0; n < m_nsurf; n++) { - m_surf[n]->setCoverages(c + loc); - loc += m_nsp[n]; - } + void ImplicitSurfChem::updateState(doublereal* c) { + int loc = 0; + for (int n = 0; n < m_nsurf; n++) { + m_surf[n]->setCoverages(c + loc); + loc += m_nsp[n]; } + } - - /** + /* * Called by the integrator to evaluate ydot given y at time 'time'. */ void ImplicitSurfChem::eval(doublereal time, doublereal* y, @@ -138,8 +200,8 @@ namespace Cantera { int loc, k, kstart; for (n = 0; n < m_nsurf; n++) { rs0 = 1.0/m_surf[n]->siteDensity(); - m_kin[n]->getNetProductionRates(DATA_PTR(m_work)); - kstart = m_kin[n]->kineticsSpeciesIndex(0,m_surfindex[n]); + m_vecKinPtrs[n]->getNetProductionRates(DATA_PTR(m_work)); + kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]); sum = 0.0; loc = 0; for (k = 1; k < m_nsp[n]; k++) { @@ -151,4 +213,182 @@ namespace Cantera { } } + // Solve for the pseudo steady-state of the surface problem + /* + * Solve for the steady state of the surface problem. + * This is the same thing as the advanceCoverages() function, + * but at infinite times. + * + * Note, a direct solve is carried out under the hood here, + * to reduce the computational time. + */ + void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride, + doublereal timeScaleOverride) { + + int ifunc; + /* + * set bulkFunc + * -> We assume that the bulk concentrations are constant. + */ + int bulkFunc = BULK_ETCH; + /* + * time scale - time over which to integrate equations + */ + double time_scale = timeScaleOverride; + /* + * + */ + if (!m_surfSolver) { + m_surfSolver = new solveSP(this, bulkFunc); + /* + * set ifunc, which sets the algorithm. + */ + ifunc = SFLUX_INITIALIZE; + } else { + ifunc = SFLUX_RESIDUAL; + } + + // Possibly override the ifunc value + if (ifuncOverride >= 0) { + ifunc = ifuncOverride; + } + + /* + * Get the specifications for the problem from the values + * in the ThermoPhase objects for all phases. + * + * 1) concentrations of all species in all phases, m_concSpecies[] + * 2) Temperature and pressure + */ + getConcSpecies(DATA_PTR(m_concSpecies)); + InterfaceKinetics *ik = m_vecKinPtrs[0]; + ThermoPhase &tp = ik->thermo(0); + double TKelvin = tp.temperature(); + double PGas = tp.pressure(); + /* + * Make sure that there is a common temperature and + * pressure for all ThermoPhase objects belonging to the + * interfacial kinetics object, if it is required by + * the problem statement. + */ + if (m_commonTempPressForPhases) { + setCommonState_TP(TKelvin, PGas); + } + + double reltol = 1.0E-6; + double atol = 1.0E-20; + + /* + * Install a filter for negative concentrations. One of the + * few ways solvess can fail is if concentrations on input + * are below zero. + */ + bool rset = false; + for (int k = 0; k < m_nv; k++) { + if (m_concSpecies[k] < 0.0) { + rset = true; + m_concSpecies[k] = 0.0; + } + } + if (rset) { + setConcSpecies(DATA_PTR(m_concSpecies)); + } + + m_surfSolver->ioflag = m_ioFlag; + + // Save the current solution + copy(m_concSpecies.begin(), m_concSpecies.end(), m_concSpeciesSave.begin()); + + + int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas, + reltol, atol); + if (retn != 1) { + // reset the concentrations + copy(m_concSpeciesSave.begin(), m_concSpeciesSave.end(), m_concSpecies.begin()); + setConcSpecies(DATA_PTR(m_concSpeciesSave)); + ifunc = SFLUX_INITIALIZE; + retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas, + reltol, atol); + + if (retn != 1) { + throw CanteraError("ImplicitSurfChem::solvePseudoSteadyStateProblem", + "solveSP return an error condition!"); + } + } + } + + + + /* + * getConcSpecies(): + * + * Fills the local concentration vector, m_concSpecies for all of the + * species in all of the phases that are unknowns in the surface + * problem. + * + * m_concSpecies[] + */ + void ImplicitSurfChem::getConcSpecies(doublereal * const vecConcSpecies) const { + int kstart; + for (int ip = 0; ip < m_nsurf; ip++) { + ThermoPhase * TP_ptr = m_surf[ip]; + kstart = m_specStartIndex[ip]; + TP_ptr->getConcentrations(vecConcSpecies + kstart); + } + kstart = m_nv; + for (int ip = 0; ip < m_numBulkPhases; ip++) { + ThermoPhase * TP_ptr = m_bulkPhases[ip]; + int nsp = TP_ptr->nSpecies(); + TP_ptr->getConcentrations(vecConcSpecies + kstart); + kstart += nsp; + } + } + + /* + * setConcSpecies(): + * + * Fills the local concentration vector, m_concSpecies for all of the + * species in all of the phases that are unknowns in the surface + * problem. + * + * m_concSpecies[] + */ + void ImplicitSurfChem::setConcSpecies(const doublereal * const vecConcSpecies) { + int kstart; + for (int ip = 0; ip < m_nsurf; ip++) { + ThermoPhase * TP_ptr = m_surf[ip]; + kstart = m_specStartIndex[ip]; + TP_ptr->setConcentrations(vecConcSpecies + kstart); + } + kstart = m_nv; + for (int ip = 0; ip < m_numBulkPhases; ip++) { + ThermoPhase * TP_ptr = m_bulkPhases[ip]; + int nsp = TP_ptr->nSpecies(); + TP_ptr->setConcentrations(vecConcSpecies + kstart); + kstart += nsp; + } + } + + /* + * setCommonState_TP(): + * + * Sets a common temperature and pressure amongst the + * thermodynamic objects in the interfacial kinetics object. + * + * Units Temperature = Kelvin + * Pressure = Pascal + */ + void ImplicitSurfChem:: + setCommonState_TP(double TKelvin, double PresPa) { + int nphases = m_nsurf; + for (int ip = 0; ip < nphases; ip++) { + ThermoPhase *TP_ptr = m_surf[ip]; + TP_ptr->setState_TP(TKelvin, PresPa); + } + for (int ip = 0; ip < m_numBulkPhases; ip++) { + ThermoPhase *TP_ptr = m_bulkPhases[ip]; + TP_ptr->setState_TP(TKelvin, PresPa); + } + } + } diff --git a/Cantera/src/kinetics/ImplicitSurfChem.h b/Cantera/src/kinetics/ImplicitSurfChem.h index bcf834071..adc0fcdce 100755 --- a/Cantera/src/kinetics/ImplicitSurfChem.h +++ b/Cantera/src/kinetics/ImplicitSurfChem.h @@ -26,9 +26,12 @@ #include "Integrator.h" #include "InterfaceKinetics.h" #include "SurfPhase.h" +#include "solveSP.h" namespace Cantera { + class solveSP; + //! Advances the surface coverages of the associated set of SurfacePhase //! objects in time @@ -115,6 +118,36 @@ namespace Cantera { * @param t1 Final Time -> This is an input */ void integrate0(doublereal t0, doublereal t1); + + + //! Solve for the pseudo steady-state of the surface problem + /*! + * Solve for the steady state of the surface problem. + * This is the same thing as the advanceCoverages() function, + * but at infinite times. + * + * Note, a direct solve is carried out under the hood here, + * to reduce the computational time. + * + * @param ifuncOverride 4 values are possible + * 1 SFLUX_INITIALIZE + * 2 SFLUX_RESIDUAL + * 3 SFLUX_JACOBIAN + * 4 SFLUX_TRANSIENT + * The default is -1, which means that the program + * will decide. + * @param timeScaleOverride When a psuedo transient is + * selected this value can be used to override + * the default time scale for integration which + * is one. + * When SFLUX_TRANSIENT is used, this is equal to the + * time over which the equations are integrated. + * When SFLUX_INITIALIZE is used, this is equal to the + * time used in the initial transient algorithm, + * before the equation system is solved directly. + */ + void solvePseudoSteadyStateProblem(int ifuncOverride = -1, + doublereal timeScaleOverride = 1.0); // overloaded methods of class FuncEval @@ -144,6 +177,59 @@ namespace Cantera { virtual void getInitialConditions(doublereal t0, size_t leny, doublereal* y); + /* + * Get the specifications for the problem from the values + * in the ThermoPhase objects for all phases. + * + * 1) concentrations of all species in all phases, m_concSpecies[] + * 2) Temperature and pressure + * + * + * @param vecConcSpecies Vector of concentrations. The + * phase concentration vectors are contiguous + * within the object, in the same order as the + * unknown vector. + */ + void getConcSpecies(doublereal * const vecConcSpecies) const; + + //! Sets the concentrations within phases that are unknowns in + //! the surface problem + /*! + * Fills the local concentration vector for all of the + * species in all of the phases that are unknowns in the surface + * problem. + * + * @param vecConcSpecies Vector of concentrations. The + * phase concentration vectors are contiguous + * within the object, in the same order as the + * unknown vector. + */ + void setConcSpecies(const doublereal * const vecConcSpecies); + + //! Sets the state variable in all thermodynamic phases (surface and + //! surrounding bulk phases) to the input temperature and pressure + /*! + * @param TKelvin input temperature (kelvin) + * @param PresPa input pressure in pascal. + */ + void setCommonState_TP(double TKelvin, double PresPa); + + + //! Returns a reference to the vector of pointers to the + //! InterfaceKinetics objects + /*! + * This should probably go away in the future, as it opens up the + * class. + */ + std::vector & getObjects() { + return m_vecKinPtrs; + } + + int checkMatch(std::vector m_vec, ThermoPhase *thPtr); + + void setIOFlag(int ioFlag) { + m_ioFlag = ioFlag; + } protected: @@ -161,20 +247,100 @@ namespace Cantera { */ void updateState(doublereal* y); + //! vector of pointers to surface phases. std::vector m_surf; - std::vector m_kin; + + //! Vector of pointers to bulk phases + std::vector m_bulkPhases; + + //! vector of pointers to InterfaceKinetics objects + std::vector m_vecKinPtrs; + + //! Vector of number of species in each Surface Phase vector_int m_nsp; + + //! index of the surface phase in each InterfaceKinetics object vector_int m_surfindex; - int m_nsurf; - int m_nv; - //int m_nsp, m_surfindex; - Integrator* m_integ; // pointer to integrator + + + vector_int m_specStartIndex; + + //! Total number of surface phases. + /*! + * This is also equal to the number of InterfaceKinetics objects + * as there is a 1-1 correspondence between InterfaceKinetics objects + * and surface phases. + */ + int m_nsurf; + + //! Total number of surface species in all surface phases + /*! + * This is the total number of unknowns in m_mode 0 problem + */ + int m_nv; + + int m_numBulkPhases; + vector_int m_nspBulkPhases; + int m_numTotalBulkSpecies; + int m_numTotalSpecies; + + std::vector pLocVec; + //! Pointer to the cvode integrator + Integrator* m_integ; doublereal m_atol, m_rtol; // tolerances doublereal m_maxstep; // max step size vector_fp m_work; + + + /** + * Temporary vector - length num species in the Kinetics object. + * This is the sum of the number of species + * in each phase included in the kinetics object. + */ + vector_fp m_concSpecies; + vector_fp m_concSpeciesSave; + + //std::vector m_vectorConcKinSpecies; + //std::vector m_vectorNetSpeciesProdRate; + /** + * Index into the species vector of the kinetics manager, + * pointing to the first species from the surrounding medium. + */ + int m_mediumSpeciesStart; + /** + * Index into the species vector of the kinetics manager, + * pointing to the first species from the condensed phase + * of the particles. + */ + int m_bulkSpeciesStart; + /** + * Index into the species vector of the kinetics manager, + * pointing to the first species from the surface + * of the particles + */ + int m_surfSpeciesStart; + /** + * Pointer to the helper method, Placid, which solves the + * surface problem. + */ + solveSP *m_surfSolver; + + //! If true, a common temperature and pressure for all + //! surface and bulk phases associated with the surface problem + //! is imposed + bool m_commonTempPressForPhases; + + //! We make the solveSS class a friend because we need + //! to access all of the above information directly. + //! Adding the members into the class is also a possibility. + friend class solveSS; + private: + //! Controls the amount of printing from this routine + //! and underlying routines. + int m_ioFlag; }; } diff --git a/Cantera/src/kinetics/InterfaceKinetics.cpp b/Cantera/src/kinetics/InterfaceKinetics.cpp index 0631d0e36..80f1c353f 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.cpp +++ b/Cantera/src/kinetics/InterfaceKinetics.cpp @@ -46,7 +46,8 @@ namespace Cantera { m_integrator(0), m_finalized(false), m_has_coverage_dependence(false), - m_has_electrochem_rxns(false) + m_has_electrochem_rxns(false), + m_ioFlag(0) { if (thermo != 0) addPhase(*thermo); m_kdata = new InterfaceKineticsData; @@ -59,7 +60,9 @@ namespace Cantera { InterfaceKinetics:: ~InterfaceKinetics(){ delete m_kdata; - delete m_integrator; + if (m_integrator) { + delete m_integrator; + } } @@ -698,6 +701,13 @@ namespace Cantera { } + void InterfaceKinetics::setIOFlag(int ioFlag) { + m_ioFlag = ioFlag; + if (m_integrator) { + m_integrator->setIOFlag(ioFlag); + } + } + // void InterfaceKinetics:: // addGlobalReaction(const ReactionData& r) { @@ -905,31 +915,39 @@ namespace Cantera { * * Note, a direct solve is carried out under the hood here, * to reduce the computational time. + * + * the integrator object is saved inbetween calls to + * reduce the computational cost of repeated calls. */ - void InterfaceKinetics::solvePseudoSteadyStateProblem() { -#ifndef DEBUG_HKM - advanceCoverages(1000.0); -#else + void InterfaceKinetics:: + solvePseudoSteadyStateProblem(int ifuncOverride, + doublereal timeScaleOverride) { + // create our own solver object + if (m_integrator == 0) { + vector k; + k.push_back(this); + m_integrator = new ImplicitSurfChem(k); + m_integrator->initialize(); + } + m_integrator->setIOFlag(m_ioFlag); /* * New direct method to go here */ - -#endif + m_integrator->solvePseudoSteadyStateProblem(ifuncOverride, timeScaleOverride); } - - void EdgeKinetics::finalize() { - m_rwork.resize(nReactions()); - int ks = reactionPhaseIndex(); - if (ks < 0) throw CanteraError("EdgeKinetics::finalize", - "no edge phase is present."); - m_surf = (SurfPhase*)&thermo(ks); - if (m_surf->nDim() != 1) - throw CanteraError("EdgeKinetics::finalize", - "expected interface dimension = 1, but got dimension = " - +int2str(m_surf->nDim())); - m_finalized = true; - } + void EdgeKinetics::finalize() { + m_rwork.resize(nReactions()); + int ks = reactionPhaseIndex(); + if (ks < 0) throw CanteraError("EdgeKinetics::finalize", + "no edge phase is present."); + m_surf = (SurfPhase*)&thermo(ks); + if (m_surf->nDim() != 1) + throw CanteraError("EdgeKinetics::finalize", + "expected interface dimension = 1, but got dimension = " + +int2str(m_surf->nDim())); + m_finalized = true; + } } diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 5aa19da09..e6bc65146 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -1,6 +1,8 @@ /** * @file InterfaceKinetics.h * + */ +/* * $Author$ * $Revision$ * $Date$ @@ -60,519 +62,546 @@ namespace Cantera { }; + /// + /// A kinetics manager for heterogeneous reaction mechanisms. The + /// reactions are assumed to occur at a 2D interface between two + /// 3D phases. + /// + class InterfaceKinetics : public Kinetics { + + public: + + /** + * Constructor + * + * @param thermo The optional parameter may be used to initialize + * the object with one ThermoPhase object. + * HKM Note -> Since the interface kinetics + * object will probably require multiple thermophase + * objects, this is probably not a good idea + * to have this parameter. + */ + InterfaceKinetics(thermo_t* thermo = 0); + + + /// Destructor. + virtual ~InterfaceKinetics(); + + virtual int ID() { return cInterfaceKinetics; } + virtual int type() { return cInterfaceKinetics; } + + /** + * Set the electric potential in the nth phase + * + * @param n phase Index in this kinetics object. + * @param V Electric potential (volts) + */ + void setElectricPotential(int n, doublereal V) { + thermo(n).setElectricPotential(V); + m_redo_rates = true; + } + + /// - /// A kinetics manager for heterogeneous reaction mechanisms. The - /// reactions are assumed to occur at a 2D interface between two - /// 3D phases. + /// @name Reaction Rates Of Progress /// - class InterfaceKinetics : public Kinetics { + //@{ - public: + //! Return the forward rates of progress for each reaction + /*! + * @param fwdROP vector of rates of progress. + * length = number of reactions, Units are kmol m-2 s-1. + */ + virtual void getFwdRatesOfProgress(doublereal* fwdROP) { + updateROP(); + std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP); + } - /** - * Constructor - * - * @param thermo The optional parameter may be used to initialize - * the object with one ThermoPhase object. - * HKM Note -> Since the interface kinetics - * object will probably require multiple thermophase - * objects, this is probably not a good idea - * to have this parameter. - */ - InterfaceKinetics(thermo_t* thermo = 0); + //! Return the reverse rates of progress for each reaction + /*! + * @param revROP vector of rates of progress. + * length = number of reactions, Units are kmol m-2 s-1. + */ + virtual void getRevRatesOfProgress(doublereal* revROP) { + updateROP(); + std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP); + } + + //! Return the net rates of progress for each reaction + /*! + * @param netROP vector of rates of progress. + * length = number of reactions, Units are kmol m-2 s-1. + */ + virtual void getNetRatesOfProgress(doublereal* netROP) { + updateROP(); + std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); + } + + virtual void getEquilibriumConstants(doublereal* kc); - /// Destructor. - virtual ~InterfaceKinetics(); + virtual void getDeltaGibbs( doublereal* deltaG); - virtual int ID() { return cInterfaceKinetics; } - virtual int type() { return cInterfaceKinetics; } - - /** - * Set the electric potential in the nth phase - * - * @param n phase Index in this kinetics object. - * @param V Electric potential (volts) - */ - void setElectricPotential(int n, doublereal V) { - thermo(n).setElectricPotential(V); - m_redo_rates = true; - } - - - /// - /// @name Reaction Rates Of Progress - /// - //@{ - - //! Return the forward rates of progress for each reaction - /*! - * @param fwdROP vector of rates of progress. - * length = number of reactions, Units are kmol m-2 s-1. - */ - virtual void getFwdRatesOfProgress(doublereal* fwdROP) { - updateROP(); - std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP); - } - - //! Return the reverse rates of progress for each reaction - /*! - * @param revROP vector of rates of progress. - * length = number of reactions, Units are kmol m-2 s-1. - */ - virtual void getRevRatesOfProgress(doublereal* revROP) { - updateROP(); - std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP); - } - - //! Return the net rates of progress for each reaction - /*! - * @param netROP vector of rates of progress. - * length = number of reactions, Units are kmol m-2 s-1. - */ - virtual void getNetRatesOfProgress(doublereal* netROP) { - updateROP(); - std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP); - } - - virtual void getEquilibriumConstants(doublereal* kc); - - - virtual void getDeltaGibbs( doublereal* deltaG); - - /** - * Return the vector of values for the reactions change in - * enthalpy. - * These values depend upon the concentration - * of the solution. - * - * units = J kmol-1 - */ - virtual void getDeltaEnthalpy( doublereal* deltaH); + /** + * Return the vector of values for the reactions change in + * enthalpy. + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 + */ + virtual void getDeltaEnthalpy( doublereal* deltaH); - //! Return the vector of values for the change in - //! entropy due to each reaction - /*! - * These values depend upon the concentration - * of the solution. - * - * units = J kmol-1 Kelvin-1 - * - * @param deltaS vector of Enthalpy changes - * Length = m_ii, number of reactions - * - */ - virtual void getDeltaEntropy(doublereal* deltaS); + //! Return the vector of values for the change in + //! entropy due to each reaction + /*! + * These values depend upon the concentration + * of the solution. + * + * units = J kmol-1 Kelvin-1 + * + * @param deltaS vector of Enthalpy changes + * Length = m_ii, number of reactions + * + */ + virtual void getDeltaEntropy(doublereal* deltaS); - //! Return the vector of values for the reaction - //! standard state gibbs free energy change. - /*! - * These values don't depend upon the concentration - * of the solution. - * - * @param deltaG vector of rxn SS free energy changes - * units = J kmol-1 - */ - virtual void getDeltaSSGibbs(doublereal* deltaG); + //! Return the vector of values for the reaction + //! standard state gibbs free energy change. + /*! + * These values don't depend upon the concentration + * of the solution. + * + * @param deltaG vector of rxn SS free energy changes + * units = J kmol-1 + */ + virtual void getDeltaSSGibbs(doublereal* deltaG); - //! Return the vector of values for the change in the - //! standard state enthalpies of reaction. - /*! - * These values don't depend upon the concentration - * of the solution. - * - * @param deltaH vector of rxn SS enthalpy changes - * units = J kmol-1 - */ - virtual void getDeltaSSEnthalpy(doublereal* deltaH); + //! Return the vector of values for the change in the + //! standard state enthalpies of reaction. + /*! + * These values don't depend upon the concentration + * of the solution. + * + * @param deltaH vector of rxn SS enthalpy changes + * units = J kmol-1 + */ + virtual void getDeltaSSEnthalpy(doublereal* deltaH); - //! Return the vector of values for the change in the - //! standard state entropies for each reaction. - /*! - * These values don't depend upon the concentration - * of the solution. - * - * @param deltaS vector of rxn SS entropy changes - * units = J kmol-1 Kelvin-1 - */ - virtual void getDeltaSSEntropy(doublereal* deltaS); + //! Return the vector of values for the change in the + //! standard state entropies for each reaction. + /*! + * These values don't depend upon the concentration + * of the solution. + * + * @param deltaS vector of rxn SS entropy changes + * units = J kmol-1 Kelvin-1 + */ + virtual void getDeltaSSEntropy(doublereal* deltaS); - //@} - /** - * @name Species Production Rates - */ - //@{ + //@} + /** + * @name Species Production Rates + */ + //@{ - //! Returns the Species creation rates [kmol/m^2/s]. - /*! - * Return the species - * creation rates in array cdot, which must be - * dimensioned at least as large as the total number of - * species in all phases of the kinetics - * model - * - * @param cdot Vector containing creation rates. - * length = m_kk. units = kmol/m^2/s - */ - virtual void getCreationRates(doublereal* cdot); + //! Returns the Species creation rates [kmol/m^2/s]. + /*! + * Return the species + * creation rates in array cdot, which must be + * dimensioned at least as large as the total number of + * species in all phases of the kinetics + * model + * + * @param cdot Vector containing creation rates. + * length = m_kk. units = kmol/m^2/s + */ + virtual void getCreationRates(doublereal* cdot); - //! Return the Species destruction rates [kmol/m^2/s]. - /*! - * Return the species destruction rates in array ddot, which must be - * dimensioned at least as large as the total number of - * species in all phases of the kinetics model - * - * @param ddot Vector containing destruction rates. - * length = m_kk. units = kmol/m^2/s - */ - virtual void getDestructionRates(doublereal* ddot); + //! Return the Species destruction rates [kmol/m^2/s]. + /*! + * Return the species destruction rates in array ddot, which must be + * dimensioned at least as large as the total number of + * species in all phases of the kinetics model + * + * @param ddot Vector containing destruction rates. + * length = m_kk. units = kmol/m^2/s + */ + virtual void getDestructionRates(doublereal* ddot); - //! Return the species net production rates [kmol/m^2/s]. - /*! - * Species net production rates [kmol/m^2/s]. Return the species - * net production rates (creation - destruction) in array - * wdot, which must be dimensioned at least as large as the - * total number of species in all phases of the kinetics - * model - * - * @param net Vector of species production rates. - * units kmol m-d s-1, where d is dimension. - */ - virtual void getNetProductionRates(doublereal* net); + //! Return the species net production rates [kmol/m^2/s]. + /*! + * Species net production rates [kmol/m^2/s]. Return the species + * net production rates (creation - destruction) in array + * wdot, which must be dimensioned at least as large as the + * total number of species in all phases of the kinetics + * model + * + * @param net Vector of species production rates. + * units kmol m-d s-1, where d is dimension. + */ + virtual void getNetProductionRates(doublereal* net); - //@} - /** - * @name Reaction Mechanism Informational Query Routines - */ - //@{ + //@} + /** + * @name Reaction Mechanism Informational Query Routines + */ + //@{ - /** - * Stoichiometric coefficient of species k as a reactant in - * reaction i. - */ - virtual doublereal reactantStoichCoeff(int k, int i) const { - return m_rrxn[k][i]; - } + /** + * Stoichiometric coefficient of species k as a reactant in + * reaction i. + */ + virtual doublereal reactantStoichCoeff(int k, int i) const { + return m_rrxn[k][i]; + } - /** - * Stoichiometric coefficient of species k as a product in - * reaction i. - */ - virtual doublereal productStoichCoeff(int k, int i) const { - return m_prxn[k][i]; - } + /** + * Stoichiometric coefficient of species k as a product in + * reaction i. + */ + virtual doublereal productStoichCoeff(int k, int i) const { + return m_prxn[k][i]; + } - /** - * Flag specifying the type of reaction. The legal values and - * their meaning are specific to the particular kinetics - * manager. - */ - virtual int reactionType(int i) const { - return m_index[i].first; - } + /** + * Flag specifying the type of reaction. The legal values and + * their meaning are specific to the particular kinetics + * manager. + */ + virtual int reactionType(int i) const { + return m_index[i].first; + } - //! Return the charge transfer rxn Beta parameter for the ith reaction - /*! - * Returns the beta parameter for a charge transfer reaction. This - * parameter is not important for non-charge transfer reactions. - * Note, the parameter defaults to zero. However, a value of 0.5 - * should be supplied for every charge transfer reaction if - * no information is known, as a value of 0.5 pertains to a - * symmetric transition state. The value can vary between 0 to 1. - * - * - * @param irxn Reaction number in the kinetics mechanism - * - * @return - * Beta parameter. This defaults to zero, even for charge transfer - * reactions. - */ - doublereal electrochem_beta(int irxn) const; + //! Return the charge transfer rxn Beta parameter for the ith reaction + /*! + * Returns the beta parameter for a charge transfer reaction. This + * parameter is not important for non-charge transfer reactions. + * Note, the parameter defaults to zero. However, a value of 0.5 + * should be supplied for every charge transfer reaction if + * no information is known, as a value of 0.5 pertains to a + * symmetric transition state. The value can vary between 0 to 1. + * + * + * @param irxn Reaction number in the kinetics mechanism + * + * @return + * Beta parameter. This defaults to zero, even for charge transfer + * reactions. + */ + doublereal electrochem_beta(int irxn) const; - /** - * True if reaction i has been declared to be reversible. If - * isReversible(i) is false, then the reverse rate of progress - * for reaction i is always zero. - */ - virtual bool isReversible(int i) { - if (std::find(m_revindex.begin(), m_revindex.end(), i) - < m_revindex.end()) return true; - else return false; - } + /** + * True if reaction i has been declared to be reversible. If + * isReversible(i) is false, then the reverse rate of progress + * for reaction i is always zero. + */ + virtual bool isReversible(int i) { + if (std::find(m_revindex.begin(), m_revindex.end(), i) + < m_revindex.end()) return true; + else return false; + } - /** - * Return a string representing the reaction. - */ - virtual std::string reactionString(int i) const { - return m_rxneqn[i]; - } + /** + * Return a string representing the reaction. + */ + virtual std::string reactionString(int i) const { + return m_rxneqn[i]; + } - virtual void getFwdRateConstants(doublereal* kfwd); - virtual void getRevRateConstants(doublereal* krev, - bool doIrreversible = false); + virtual void getFwdRateConstants(doublereal* kfwd); + virtual void getRevRateConstants(doublereal* krev, + bool doIrreversible = false); - virtual void getActivationEnergies(doublereal *E); + virtual void getActivationEnergies(doublereal *E); - //@} - /** - * @name Reaction Mechanism Construction - */ - //@{ + //@} + /** + * @name Reaction Mechanism Construction + */ + //@{ - //! Prepare the class for the addition of reactions. - /*! - * This function must be called after instantiation of the class, but before - * any reactions are actually added to the mechanism. - * This function calculates m_kk the number of species in all - * phases participating in the reaction mechanism. We don't know - * m_kk previously, before all phases have been added. - */ - virtual void init(); + //! Prepare the class for the addition of reactions. + /*! + * This function must be called after instantiation of the class, but before + * any reactions are actually added to the mechanism. + * This function calculates m_kk the number of species in all + * phases participating in the reaction mechanism. We don't know + * m_kk previously, before all phases have been added. + */ + virtual void init(); - //! Add a single reaction to the mechanism. - /*! - * @param r Reference to a ReactionData object containing all of - * the info needed to describe the reaction. - */ - virtual void addReaction(const ReactionData& r); + //! Add a single reaction to the mechanism. + /*! + * @param r Reference to a ReactionData object containing all of + * the info needed to describe the reaction. + */ + virtual void addReaction(const ReactionData& r); - //! Finish adding reactions and prepare for use. - /*! - * This function - * must be called after all reactions are entered into the mechanism - * and before the mechanism is used to calculate reaction rates. - */ - virtual void finalize(); + //! Finish adding reactions and prepare for use. + /*! + * This function + * must be called after all reactions are entered into the mechanism + * and before the mechanism is used to calculate reaction rates. + */ + virtual void finalize(); - virtual bool ready() const; + virtual bool ready() const; - void updateROP(); + void updateROP(); - void _update_rates_T(); - void _update_rates_phi(); - void _update_rates_C(); + void _update_rates_T(); + void _update_rates_phi(); + void _update_rates_C(); - //! Advance the surface coverages in time - /*! - * This method carries out a time-accurate advancement of the - * surface coverages for a specified amount of time. - * - * \f[ - * \dot {\theta}_k = \dot s_k (\sigma_k / s_0) - * \f] - * - * - * @param tstep Time value to advance the surface coverages - */ - void advanceCoverages(doublereal tstep); + //! Advance the surface coverages in time + /*! + * This method carries out a time-accurate advancement of the + * surface coverages for a specified amount of time. + * + * \f[ + * \dot {\theta}_k = \dot s_k (\sigma_k / s_0) + * \f] + * + * + * @param tstep Time value to advance the surface coverages + */ + void advanceCoverages(doublereal tstep); - //! Solve for the pseudo steady-state of the surface problem - /*! - * Solve for the steady state of the surface problem. - * This is the same thing as the advanceCoverages() function, - * but at infinite times. - * - * Note, a direct solve is carried out under the hood here, - * to reduce the computational time. - */ - void solvePseudoSteadyStateProblem(); + //! Solve for the pseudo steady-state of the surface problem + /*! + * Solve for the steady state of the surface problem. + * This is the same thing as the advanceCoverages() function, + * but at infinite times. + * + * Note, a direct solve is carried out under the hood here, + * to reduce the computational time. + * + * @param ifuncOverride 4 values are possible + * 1 SFLUX_INITIALIZE + * 2 SFLUX_RESIDUAL + * 3 SFLUX_JACOBIAN + * 4 SFLUX_TRANSIENT + * The default is -1, which means that the program + * will decide. + * @param timeScaleOverride When a psuedo transient is + * selected this value can be used to override + * the default time scale for integration which + * is one. + * When SFLUX_TRANSIENT is used, this is equal to the + * time over which the equations are integrated. + * When SFLUX_INITIALIZE is used, this is equal to the + * time used in the initial transient algorithm, + * before the equation system is solved directly. + */ + void solvePseudoSteadyStateProblem(int ifuncOverride = -1, + doublereal timeScaleOverride = 1.0); - void checkPartialEquil(); + void setIOFlag(int ioFlag); + + void checkPartialEquil(); - //! Temporary work vector of length m_kk - vector_fp m_grt; + //! Temporary work vector of length m_kk + vector_fp m_grt; - protected: + protected: - //! m_kk is the number of species in all of the phases - //! that participate in this kinetics mechanism. - int m_kk; + //! m_kk is the number of species in all of the phases + //! that participate in this kinetics mechanism. + int m_kk; - //! List of reactions numbers which are reversible reactions - /*! - * This is a vector of reaction numbers. Each reaction - * in the list is reversible. - * Length = number of reversible reactions - */ - vector_int m_revindex; + //! List of reactions numbers which are reversible reactions + /*! + * This is a vector of reaction numbers. Each reaction + * in the list is reversible. + * Length = number of reversible reactions + */ + vector_int m_revindex; - Rate1 m_rates; - bool m_redo_rates; + Rate1 m_rates; + bool m_redo_rates; - /** - * Vector of information about reactions in the - * mechanism. - * The key is the reaction index (0 < i < m_ii). - * The first pair is the reactionType of the reaction. - * The second pair is ... - */ - mutable std::map > m_index; + /** + * Vector of information about reactions in the + * mechanism. + * The key is the reaction index (0 < i < m_ii). + * The first pair is the reactionType of the reaction. + * The second pair is ... + */ + mutable std::map > m_index; - //! Vector of irreversible reaction numbers - /*! - * vector containing the reaction numbers of irreversible - * reactions. - */ - std::vector m_irrev; + //! Vector of irreversible reaction numbers + /*! + * vector containing the reaction numbers of irreversible + * reactions. + */ + std::vector m_irrev; - //! Stoichiometric manager for the reaction mechanism - /*! - * This is the manager for the kinetics mechanism that - * handles turning reaction extents into species - * production rates and also handles turning thermo - * properties into reaction thermo properties. - */ - ReactionStoichMgr m_rxnstoich; + //! Stoichiometric manager for the reaction mechanism + /*! + * This is the manager for the kinetics mechanism that + * handles turning reaction extents into species + * production rates and also handles turning thermo + * properties into reaction thermo properties. + */ + ReactionStoichMgr m_rxnstoich; - //! Number of irreversible reactions in the mechanism - int m_nirrev; + //! Number of irreversible reactions in the mechanism + int m_nirrev; - //! Number of reversible reactions in the mechanism - int m_nrev; + //! Number of reversible reactions in the mechanism + int m_nrev; - //! m_rrxn is a vector of maps, containing the reactant - //! stochiometric coefficient information - /*! - * m_rrxn has a length - * equal to the total number of species in the kinetics - * object. For each species, there exists a map, with the - * reaction number being the key, and the - * reactant stoichiometric coefficient for the species being the value. - * HKM -> mutable because search sometimes creates extra - * entries. To be fixed in future... - */ - mutable std::vector > m_rrxn; + //! m_rrxn is a vector of maps, containing the reactant + //! stochiometric coefficient information + /*! + * m_rrxn has a length + * equal to the total number of species in the kinetics + * object. For each species, there exists a map, with the + * reaction number being the key, and the + * reactant stoichiometric coefficient for the species being the value. + * HKM -> mutable because search sometimes creates extra + * entries. To be fixed in future... + */ + mutable std::vector > m_rrxn; - //! m_prxn is a vector of maps, containing the reactant - //! stochiometric coefficient information - /** - * m_prxn is a vector of maps. m_prxn has a length - * equal to the total number of species in the kinetics - * object. For each species, there exists a map, with the - * reaction number being the key, and the - * product stoichiometric coefficient for the species being the value. - */ - mutable std::vector > m_prxn; - //! String expression for each rxn - /*! - * Vector of strings of length m_ii, the number of - * reactions, containing the - * string expressions for each reaction - * (e.g., reactants <=> product1 + product2) - */ - std::vector m_rxneqn; + //! m_prxn is a vector of maps, containing the reactant + //! stochiometric coefficient information + /** + * m_prxn is a vector of maps. m_prxn has a length + * equal to the total number of species in the kinetics + * object. For each species, there exists a map, with the + * reaction number being the key, and the + * product stoichiometric coefficient for the species being the value. + */ + mutable std::vector > m_prxn; + //! String expression for each rxn + /*! + * Vector of strings of length m_ii, the number of + * reactions, containing the + * string expressions for each reaction + * (e.g., reactants <=> product1 + product2) + */ + std::vector m_rxneqn; - /** - * Temporary data storage used in calculating the rates of - * of reactions. - */ - InterfaceKineticsData* m_kdata; + /** + * Temporary data storage used in calculating the rates of + * of reactions. + */ + InterfaceKineticsData* m_kdata; - //! an array of generalized concentrations for each species - /*! - * An array of generalized concentrations - * \f$ C_k \f$ that are defined such that \f$ a_k = C_k / - * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/ - * These generalized concentrations are used - * by this kinetics manager class to compute the forward and - * reverse rates of elementary reactions. The "units" for the - * concentrations of each phase depend upon the implementation - * of kinetics within that phase. - * The order of the species within the vector is based on - * the order of listed ThermoPhase objects in the class, and the - * order of the species within each ThermoPhase class. - */ - vector_fp m_conc; + //! an array of generalized concentrations for each species + /*! + * An array of generalized concentrations + * \f$ C_k \f$ that are defined such that \f$ a_k = C_k / + * C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration/ + * These generalized concentrations are used + * by this kinetics manager class to compute the forward and + * reverse rates of elementary reactions. The "units" for the + * concentrations of each phase depend upon the implementation + * of kinetics within that phase. + * The order of the species within the vector is based on + * the order of listed ThermoPhase objects in the class, and the + * order of the species within each ThermoPhase class. + */ + vector_fp m_conc; - //! Vector of standard state chemical potentials - /*! - * This vector contains a temporary vector of - * standard state chemical potentials - * for all of the species in the kinetics object - * - * Length = m_k - * units = J/kmol - */ - vector_fp m_mu0; + //! Vector of standard state chemical potentials + /*! + * This vector contains a temporary vector of + * standard state chemical potentials + * for all of the species in the kinetics object + * + * Length = m_k + * units = J/kmol + */ + vector_fp m_mu0; - //! Vector of phase potentials - /*! - * Temporary vector containing the potential of each phase - * in the kinetics object - * - * length = number of phases - * units = Volts - */ - vector_fp m_phi; + //! Vector of phase potentials + /*! + * Temporary vector containing the potential of each phase + * in the kinetics object + * + * length = number of phases + * units = Volts + */ + vector_fp m_phi; - //! Vector of potential energies due to Voltages - /*! - * Length is the number of species in kinetics mech. It's - * used to store the potential energy due to the voltage. - */ - vector_fp m_pot; + //! Vector of potential energies due to Voltages + /*! + * Length is the number of species in kinetics mech. It's + * used to store the potential energy due to the voltage. + */ + vector_fp m_pot; - //! Vector temporary - /*! - * Length is number of reactions. it's used to store the - * voltage contribution to the activation energy. - */ - vector_fp m_rwork; + //! Vector temporary + /*! + * Length is number of reactions. it's used to store the + * voltage contribution to the activation energy. + */ + vector_fp m_rwork; - //! Vector of raw activation energies for the reactions - /*! - * units are in Kelvin - */ - vector_fp m_E; + //! Vector of raw activation energies for the reactions + /*! + * units are in Kelvin + */ + vector_fp m_E; - //! Pointer to the surface phase - SurfPhase* m_surf; + //! Pointer to the single surface phase + SurfPhase* m_surf; - //! Pointer to the surface solver - ImplicitSurfChem* m_integrator; + //! Pointer to the Implicit surface chemistry object + /*! + * Note this object is owned by this InterfaceKinetics + * object. It may only be used to solve this single + * InterfaceKinetics objects's surface problem uncoupled + * from other surface phases. + */ + ImplicitSurfChem* m_integrator; - vector_fp m_beta; - vector_int m_ctrxn; + vector_fp m_beta; + vector_int m_ctrxn; - int reactionNumber(){ return m_ii;} + int reactionNumber(){ return m_ii;} - void addElementaryReaction(const ReactionData& r); - void addGlobalReaction(const ReactionData& r); - void installReagents(const ReactionData& r); + void addElementaryReaction(const ReactionData& r); + void addGlobalReaction(const ReactionData& r); + void installReagents(const ReactionData& r); - void updateKc(); + void updateKc(); - //! Write values into m_index - /*! - * @param rxnNumber reaction number - * @param type reaction type - * @param loc location ?? - */ - void registerReaction(int rxnNumber, int type, int loc) { - m_index[rxnNumber] = std::pair(type, loc); - } + //! Write values into m_index + /*! + * @param rxnNumber reaction number + * @param type reaction type + * @param loc location ?? + */ + void registerReaction(int rxnNumber, int type, int loc) { + m_index[rxnNumber] = std::pair(type, loc); + } - void applyButlerVolmerCorrection(doublereal* kf); + void applyButlerVolmerCorrection(doublereal* kf); - //! boolean indicating whether mechanism has been finalized - bool m_finalized; - bool m_has_coverage_dependence; - bool m_has_electrochem_rxns; + //! boolean indicating whether mechanism has been finalized + bool m_finalized; + bool m_has_coverage_dependence; + bool m_has_electrochem_rxns; - private: + int m_ioFlag; + private: - }; + }; } #endif diff --git a/Cantera/src/kinetics/Makefile.in b/Cantera/src/kinetics/Makefile.in index 047cefafb..3e4ad8f6b 100644 --- a/Cantera/src/kinetics/Makefile.in +++ b/Cantera/src/kinetics/Makefile.in @@ -27,7 +27,7 @@ else endif PIC_FLAG=@PIC@ - +LOCAL_DEFS = -DDEBUG_SOLVESP CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) # homogeneous kinetics @@ -35,12 +35,12 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG) ifeq ($(do_kinetics),1) KINETICS_OBJ=importKinetics.o GRI_30_Kinetics.o KineticsFactory.o \ GasKinetics.o \ - FalloffFactory.o ReactionStoichMgr.o Kinetics.o + FalloffFactory.o ReactionStoichMgr.o Kinetics.o solveSP.o KINETICS_H = importKinetics.h GRI_30_Kinetics.h KineticsFactory.h \ Kinetics.h GasKinetics.h \ FalloffFactory.h ReactionStoichMgr.h reaction_defs.h \ FalloffMgr.h ThirdBodyMgr.h RateCoeffMgr.h ReactionData.h \ - RxnRates.h Enhanced3BConc.h StoichManager.h + RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h KINETICS = $(KINETICS_OBJ) $(KINETICS_H) endif diff --git a/Cantera/src/kinetics/solveSP.cpp b/Cantera/src/kinetics/solveSP.cpp new file mode 100644 index 000000000..dd929f271 --- /dev/null +++ b/Cantera/src/kinetics/solveSP.cpp @@ -0,0 +1,1445 @@ +/* + * @file: solveSP.cpp Implicit surface site concentration solver + */ +/* + * $Id$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#include "solveSP.h" +#include "clockWC.h" +#include "ctlapack.h" + +/* Standard include files */ + +#include +#include +#include + +#include + +using namespace std; +namespace Cantera { + + /*************************************************************************** + * STATIC ROUTINES DEFINED IN THIS FILE + ***************************************************************************/ + + static double calc_damping(double *x, double *dx, int dim, int *); + static double calcWeightedNorm(const double [], const double dx[], int); + + /*************************************************************************** + * LAPACK PROTOTYPES + ***************************************************************************/ +#define FSUB_TYPE void + extern "C" { + extern FSUB_TYPE dgetrf_(int *, int *, double *, int *, int [], int *); + extern FSUB_TYPE dgetrs_(char *, int *, int *, double *, int *, int [], + double [], int *, int *, unsigned int); + } + /***************************************************************************** + * PROTOTYPES and PREPROC DIRECTIVES FOR MISC. ROUTINES + *****************************************************************************/ + +#ifndef MAX +# define MAX(x,y) (( (x) > (y) ) ? (x) : (y)) /* max function */ +#endif + +#ifndef DAMPING +# define DAMPING true +#endif + + /*************************************************************************** + * solveSP Class Definitinos + ***************************************************************************/ + + // Main constructor + solveSP::solveSP(ImplicitSurfChem *surfChemPtr, int bulkFunc) : + m_SurfChemPtr(surfChemPtr), + m_objects(surfChemPtr->getObjects()), + m_neq(0), + m_bulkFunc(bulkFunc), + m_numSurfPhases(0), + m_numTotSurfSpecies(0), + m_numBulkPhasesSS(0), + m_numTotBulkSpeciesSS(0), + m_atol(1.0E-15), + m_rtol(1.0E-4), + m_maxstep(1000), + m_maxTotSpecies(0), + ioflag(0) + { + + m_numSurfPhases = 0; + int numPossibleSurfPhases = m_objects.size(); + for (int n = 0; n < numPossibleSurfPhases; n++) { + InterfaceKinetics *m_kin = m_objects[n]; + int surfPhaseIndex = m_kin->surfacePhaseIndex(); + if (surfPhaseIndex >= 0) { + m_numSurfPhases++; + m_indexKinObjSurfPhase.push_back(n); + m_kinObjPhaseIDSurfPhase.push_back(surfPhaseIndex); + } else { + throw CanteraError("solveSP", + "InterfaceKinetics object has no surface phase"); + } + ThermoPhase *tp = &(m_kin->thermo(surfPhaseIndex)); + SurfPhase *sp = dynamic_cast(tp); + if (!sp) { + throw CanteraError("solveSP", + "Inconsistent ThermoPhase object within " + "InterfaceKinetics object"); + } + + m_ptrsSurfPhase.push_back(sp); + int nsp = sp->nSpecies(); + m_nSpeciesSurfPhase.push_back(nsp); + m_numTotSurfSpecies += nsp; + + } + /* + * We rely on ordering to figure things out + */ + if (1) { + //m_numBulkPhases = m_kin0->nPhases() - 1 - m_numSurfPhases; + // Disable the capability until we figure out what is going on + m_numBulkPhasesSS = 0; + //if (m_numBulkPhasesSS > 0) { + //m_numBulkSpecies.resize(m_numBulkPhasesSS, 0); + //m_bulkPhasePtrs.resize(m_numBulkPhasesSS, 0); + //m_bulkIndex = 1; + //if (m_bulkIndex == surfPhaseIndex) { + // m_bulkIndex += m_numSurfPhases; + //} + + //for (i = 0; i < m_numBulkPhasesSS; i++) { + // m_bulkPhasePtrs[i] = &(m_kin0->thermo(m_bulkIndex + i)); + // m_numBulkSpecies[i] = m_bulkPhasePtrs[i]->nSpecies(); + // m_numTotBulkSpeciesSS += m_numBulkSpecies[i]; + //} + //} + } + + if (bulkFunc == BULK_DEPOSITION) { + m_neq = m_numTotSurfSpecies + m_numTotBulkSpeciesSS; + } else { + m_neq = m_numTotSurfSpecies; + } + + m_maxTotSpecies = 0; + for (int n = 0; n < m_numSurfPhases; n++) { + int tsp = m_objects[n]->nTotalSpecies(); + m_maxTotSpecies = MAX(m_maxTotSpecies, tsp); + } + m_maxTotSpecies = MAX(m_maxTotSpecies, m_neq); + + + m_netProductionRatesSave.resize(m_maxTotSpecies, 0.0); + m_numEqn1.resize(m_maxTotSpecies, 0.0); + m_numEqn2.resize(m_maxTotSpecies, 0.0); + m_XMolKinSpecies.resize(m_maxTotSpecies, 0.0); + m_CSolnSave.resize(m_neq, 0.0); + + m_spSurfLarge.resize(m_numSurfPhases, 0); + + m_kinSpecIndex.resize(m_numTotSurfSpecies + m_numTotBulkSpeciesSS, 0); + m_kinObjIndex.resize(m_numTotSurfSpecies + m_numTotBulkSpeciesSS, 0); + m_eqnIndexStartSolnPhase.resize(m_numSurfPhases + m_numBulkPhasesSS, 0); + + int kindexSP = 0; + int isp, k, nsp, kstart; + for (isp = 0; isp < m_numSurfPhases; isp++) { + int iKinObject = m_indexKinObjSurfPhase[isp]; + InterfaceKinetics *m_kin = m_objects[iKinObject]; + int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; + kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex); + nsp = m_nSpeciesSurfPhase[isp]; + m_eqnIndexStartSolnPhase[isp] = kindexSP; + for (k = 0; k < nsp; k++, kindexSP++) { + m_kinSpecIndex[kindexSP] = kstart + k; + m_kinObjIndex[kindexSP] = isp; + } + } + if (0) { + //for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + //nt iKinObject = m_bulkKinObjID[isp]; + //InterfaceKinetics *m_kin = m_objects[iKinObject]; + //int bulkIndex = m_bulkKinObjPhaseID[isp]; + //kstart = m_kin->kineticsSpeciesIndex(0, bulkIndex); + // nsp = m_numBulkSpecies[isp]; + //m_eqnIndexStartSolnPhase[isp] = kindexSP; + //for (k = 0; k < nsp; k++, kindexSP++) { + // m_kinSpecIndex[kindexSP] = kstart + k; + // m_kinObjIndex[kindexSP] = m_numSurfPhases + isp; + //} + //} + } + + // Dimension solution vector + int dim1 = MAX(1, m_neq); + m_CSolnSP.resize(dim1, 0.0); + m_CSolnSPInit.resize(dim1, 0.0); + m_CSolnSPOld.resize(dim1, 0.0); + m_wtResid.resize(dim1, 0.0); + m_wtSpecies.resize(dim1, 0.0); + m_resid.resize(dim1, 0.0); + m_ipiv.resize(dim1, 0); + + m_Jac.resize(dim1, dim1, 0.0); + m_JacCol.resize(dim1, 0); + for (int k = 0; k < dim1; k++) { + m_JacCol[k] = m_Jac.ptrColumn(k); + } + } + + // Empty destructor + solveSP::~solveSP() { + } + + /* + * The following calculation is a Newton's method to + * get the surface fractions of the surface and bulk species by + * requiring that the + * surface species production rate = 0 and that the bulk fractions are + * proportional to their production rates. + */ + int solveSP::solveSurfProb(int ifunc, double time_scale, double TKelvin, + double PGas, double reltol, double abstol) + { + double EXTRA_ACCURACY = 0.001; + if (ifunc == SFLUX_JACOBIAN) { + EXTRA_ACCURACY *= 0.001; + } + int k, irow; + int jcol, info = 0; + int label_t=-1; /* Species IDs for time control */ + int label_d; /* Species IDs for damping control */ + int label_t_old=-1; + double label_factor = 1.0; + int iter=0; // iteration number on numlinear solver + int iter_max=1000; // maximum number of nonlinear iterations + int nrhs=1; + double deltaT = 1.0E-10; // Delta time step + double damp=1.0, tmp; + // Weighted L2 norm of the residual. Currently, this is only + // used for IO purposes. It doesn't control convergence. + // Therefore, it is turned off when DEBUG_SOLVESP isn't defined. + double resid_norm; + double inv_t = 0.0; + double t_real = 0.0, update_norm = 1.0E6; + + bool do_time = false, not_converged = true; + +#ifdef DEBUG_SOLVESP +#ifdef DEBUG_SOLVESP_TIME + double t1; +#endif +#else + if (ioflag > 1) { + ioflag = 1; + } +#endif + +#ifdef DEBUG_SOLVESP +#ifdef DEBUG_SOLVESP_TIME + Cantera::clockWC wc; + if (ioflag) t1 = wc.secondsWC(); +#endif +#endif + + /* + * Set the initial value of the do_time parameter + */ + if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) do_time = true; + + /* + * Store the initial guess for the surface problem in the soln vector, + * CSoln, and in an separate vector CSolnInit. + */ + int loc = 0; + for (int n = 0; n < m_numSurfPhases; n++) { + SurfPhase *sf_ptr = m_ptrsSurfPhase[n]; + sf_ptr->getConcentrations(DATA_PTR(m_numEqn1)); + int nsp = m_nSpeciesSurfPhase[n]; + for (k = 0; k getConcentrations(DATA_PTR(m_numEqn1)); + //int nsp = m_numBulkSpecies[isp]; + //for (k = 0; k < nsp; k++, kindex++) { + // m_CSolnSP[loc] = m_numEqn1[k]; + // loc++; + //} + //} + } + std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPInit.begin()); + + // Calculate the largest species in each phase + evalSurfLarge(DATA_PTR(m_CSolnSP)); + /* + * Get the net production rate of all species in the kinetics manager. + */ + // m_kin->getNetProductionRates(DATA_PTR(m_netProductionRatesSave)); + + if (ioflag) { + print_header(ioflag, ifunc, time_scale, DAMPING, reltol, abstol, + TKelvin, PGas, DATA_PTR(m_netProductionRatesSave), + DATA_PTR(m_XMolKinSpecies)); + } + + /* + * Quick return when there isn't a surface problem to solve + */ + if (m_neq == 0) { + not_converged = false; + update_norm = 0.0; + } + + /* ------------------------------------------------------------------ + * Start of Newton's method + * ------------------------------------------------------------------ + */ + while (not_converged && iter < iter_max) { + iter++; + /* + * Store previous iteration's solution in the old solution vector + */ + std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin()); + + /* + * Evaluate the largest surface species for each surface phase every + * 5 iterations. + */ + if (iter%5 == 4) { + evalSurfLarge(DATA_PTR(m_CSolnSP)); + } + + /* + * Calculate the value of the time step + * - heuristics to stop large oscillations in deltaT + */ + if (do_time) { + /* don't hurry increase in time step at the same time as damping */ + if (damp < 1.0) label_factor = 1.0; + tmp = calc_t(DATA_PTR(m_netProductionRatesSave), + DATA_PTR(m_XMolKinSpecies), + &label_t, &label_t_old, &label_factor, ioflag); + if (iter < 10) + inv_t = tmp; + else if (tmp > 2.0*inv_t) + inv_t = 2.0*inv_t; + else { + inv_t = tmp; + } + + /* + * Check end condition + */ + + if (ifunc == SFLUX_TRANSIENT) { + tmp = t_real + 1.0/inv_t; + if (tmp > time_scale) inv_t = 1.0/(time_scale - t_real); + } + } + else { + /* make steady state calc a step of 1 million seconds to + prevent singular jacobians for some pathological cases */ + inv_t = 1.0e-6; + } + deltaT = 1.0/inv_t; + + /* + * Call the routine to numerically evaluation the jacobian + * and residual for the current iteration. + */ + resjac_eval(m_JacCol, DATA_PTR(m_resid), DATA_PTR(m_CSolnSP), + DATA_PTR(m_CSolnSPOld), do_time, deltaT); + + /* + * Calculate the weights. Make sure the calculation is carried + * out on the first iteration. + */ + if (iter%4 == 1) { + calcWeights(DATA_PTR(m_wtSpecies), DATA_PTR(m_wtResid), + m_Jac, DATA_PTR(m_CSolnSP), abstol, reltol); + } + + /* + * Find the weighted norm of the residual + */ + resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid), + DATA_PTR(m_resid), m_neq); + +#ifdef DEBUG_SOLVESP + if (ioflag > 1) { + printIterationHeader(ioflag, damp, inv_t, t_real, iter, do_time); + /* + * Print out the residual and jacobian + */ + printResJac(ioflag, m_neq, m_Jac, DATA_PTR(m_resid), + DATA_PTR(m_wtResid), resid_norm); + } +#endif + + /* + * Solve Linear system (with LAPACK). The solution is in resid[] + */ + // (void) dgetrf_(&m_neq, &m_neq, m_JacCol[0], &m_neq, + // DATA_PTR(m_ipiv), &info); + ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info); + if (info==0) { + //(void) dgetrs_(&cflag, &m_neq, &nrhs, m_JacCol[0], &m_neq, + // DATA_PTR(m_ipiv), DATA_PTR(m_resid), &m_neq, + // &info, 1); + ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0], m_neq, + DATA_PTR(m_ipiv), DATA_PTR(m_resid), m_neq, + info); + } + /* + * Force convergence if residual is small to avoid + * "nan" results from the linear solve. + */ + else { + if (ioflag) { + printf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n", + resid_norm, info); + } + for (jcol = 0; jcol < m_neq; jcol++) m_resid[jcol] = 0.0; + + /* print out some helpful info */ + if (ioflag > 1) { + printf("-----\n"); + printf("solveSurfProb: iter %d t_real %g delta_t %g\n\n", + iter,t_real, 1.0/inv_t); + printf("solveSurfProb: init guess, current concentration," + "and prod rate:\n"); + for (jcol = 0; jcol < m_neq; jcol++) { + printf("\t%d %g %g %g\n", jcol, m_CSolnSPInit[jcol], m_CSolnSP[jcol], + m_netProductionRatesSave[m_kinSpecIndex[jcol]]); + } + printf("-----\n"); + } + if (do_time) t_real += time_scale; +#ifdef DEBUG_SOLVESP + if (ioflag) { + printf("\nResidual is small, forcing convergence!\n"); + } +#endif + } + + /* + * Calculate the Damping factor needed to keep all unknowns + * between 0 and 1, and not allow too large a change (factor of 2) + * in any unknown. + */ + +#ifdef DAMPING + damp = calc_damping( DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d); +#endif + + /* + * Calculate the weighted norm of the update vector + * Here, resid is the delta of the solution, in concentration + * units. + */ + update_norm = calcWeightedNorm(DATA_PTR(m_wtSpecies), + DATA_PTR(m_resid), m_neq); + /* + * Update the solution vector and real time + * Crop the concentrations to zero. + */ + for (irow = 0; irow < m_neq; irow++) m_CSolnSP[irow] -= damp * m_resid[irow]; + for (irow = 0; irow < m_neq; irow++) { + m_CSolnSP[irow] = MAX(0.0, m_CSolnSP[irow]); + } + updateState( DATA_PTR(m_CSolnSP)); + + if (do_time) t_real += damp/inv_t; + + if (ioflag) { + printIteration(ioflag, damp, label_d, label_t, inv_t, t_real, iter, + update_norm, resid_norm, + DATA_PTR(m_netProductionRatesSave), + DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), + DATA_PTR(m_XMolKinSpecies), DATA_PTR(m_wtSpecies), + m_neq, do_time); + } + + if (ifunc == SFLUX_TRANSIENT) + not_converged = (t_real < time_scale); + else { + if (do_time) { + if (t_real > time_scale || + (resid_norm < 1.0e-7 && + update_norm*time_scale/t_real < EXTRA_ACCURACY) ) { + do_time = false; +#ifdef DEBUG_SOLVESP + if (ioflag > 1) { + printf("\t\tSwitching to steady solve.\n"); + } +#endif + } + } + else { + not_converged = ((update_norm > EXTRA_ACCURACY) || + (resid_norm > EXTRA_ACCURACY)); + } + } + } /* End of Newton's Method while statement */ + + /* + * End Newton's method. If not converged, print error message and + * recalculate sdot's at equal site fractions. + */ + if (not_converged) { + if (ioflag) { + printf("#$#$#$# Error in solveSP $#$#$#$ \n"); + printf("Newton iter on surface species did not converge, " + "update_norm = %e \n", update_norm); + printf("Continuing anyway\n"); + } + } +#ifdef DEBUG_SOLVESP +#ifdef DEBUG_SOLVESP_TIME + if (ioflag) { + printf("\nEnd of solve, time used: %e\n", wc.secondsWC()-t1); + } +#endif +#endif + + /* + * Decide on what to return in the solution vector + * - right now, will always return the last solution + * no matter how bad + */ + if (ioflag) { + fun_eval(DATA_PTR(m_resid), DATA_PTR(m_CSolnSP), DATA_PTR(m_CSolnSPOld), + false, deltaT); + resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid), + DATA_PTR(m_resid), m_neq); + printFinal(ioflag, damp, label_d, label_t, inv_t, t_real, iter, + update_norm, resid_norm, DATA_PTR(m_netProductionRatesSave), + DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), + DATA_PTR(m_XMolKinSpecies), DATA_PTR(m_wtSpecies), + DATA_PTR(m_wtResid), m_neq, do_time, + TKelvin, PGas); + } + + /* + * Return with the appropriate flag + */ + if (update_norm > 1.0) { + return -1; + } + return 1; + } + +#undef DAMPING + + + /* + * Update the surface states of the surface phases. + */ + void solveSP::updateState(const double *CSolnSP) { + int loc = 0; + for (int n = 0; n < m_numSurfPhases; n++) { + m_ptrsSurfPhase[n]->setConcentrations(CSolnSP + loc); + loc += m_nSpeciesSurfPhase[n]; + } + //if (m_bulkFunc == BULK_DEPOSITION) { + // for (int n = 0; n < m_numBulkPhasesSS; n++) { + // m_bulkPhasePtrs[n]->setConcentrations(CSolnSP + loc); + // loc += m_numBulkSpecies[n]; + // } + //} + } + + /* + * Update the mole fractions for phases which are part of the equation set + */ + void solveSP::updateMFSolnSP(double *XMolSolnSP) { + for (int isp = 0; isp < m_numSurfPhases; isp++) { + int keqnStart = m_eqnIndexStartSolnPhase[isp]; + m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart); + } + //if (m_bulkFunc == BULK_DEPOSITION) { + // for (int isp = 0; isp < m_numBulkPhasesSS; isp++) { + // int keqnStart = m_eqnIndexStartSolnPhase[isp + m_numSurfPhases]; + // m_bulkPhasePtrs[isp]->getMoleFractions(XMolSolnSP + keqnStart); + // } + //} + } + + /* + * Update the mole fractions for phases which are part of a single + * interfacial kinetics object + */ + void solveSP::updateMFKinSpecies(double *XMolKinSpecies, int isp) { + InterfaceKinetics *m_kin = m_objects[isp]; + int nph = m_kin->nPhases(); + for (int iph = 0; iph < nph; iph++) { + int ksi = m_kin->kineticsSpeciesIndex(0, iph); + ThermoPhase &thref = m_kin->thermo(iph); + thref.getMoleFractions(XMolKinSpecies + ksi); + } + } + + /* + * Update the vector that keeps track of the largest species in each + * surface phase. + */ + void solveSP::evalSurfLarge(const double *CSolnSP) { + int kindexSP = 0; + for (int isp = 0; isp < m_numSurfPhases; isp++) { + int nsp = m_nSpeciesSurfPhase[isp]; + double Clarge = CSolnSP[kindexSP]; + m_spSurfLarge[isp] = 0; + kindexSP++; + for (int k = 1; k < nsp; k++, kindexSP++) { + if (CSolnSP[kindexSP] > Clarge) { + Clarge = CSolnSP[kindexSP]; + m_spSurfLarge[isp] = k; + } + } + } + } + + /* + * This calculates the net production rates of all species + * + * This calculates the function eval. + * (should switch to special_species formulation for sum condition) + * + * @internal + * This routine uses the m_numEqn1 and m_netProductionRatesSave vectors + * as temporary internal storage. + */ + void solveSP::fun_eval(double* resid, const double *CSoln, + const double *CSolnOld, const bool do_time, + const double deltaT) + { + int isp, nsp, kstart, k, kindexSP, kins, kspecial; + double lenScale = 1.0E-9; + double sd = 0.0; + double grRate; + if (m_numSurfPhases > 0) { + /* + * update the surface concentrations with the input surface + * concentration vector + */ + updateState(CSoln); + /* + * Get the net production rates of all of the species in the + * surface kinetics mechanism + * + * HKM Should do it here for all kinetics objects so that + * bulk will eventually work. + */ + + if (do_time) { + kindexSP = 0; + for (isp = 0; isp < m_numSurfPhases; isp++) { + nsp = m_nSpeciesSurfPhase[isp]; + InterfaceKinetics *kinPtr = m_objects[isp]; + int surfIndex = kinPtr->surfacePhaseIndex(); + kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); + kins = kindexSP; + kinPtr->getNetProductionRates(DATA_PTR(m_netProductionRatesSave)); + for (k = 0; k < nsp; k++, kindexSP++) { + resid[kindexSP] = + (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT + - m_netProductionRatesSave[kstart + k]; + } + + kspecial = kins + m_spSurfLarge[isp]; + sd = m_ptrsSurfPhase[isp]->siteDensity(); + resid[kspecial] = sd; + for (k = 0; k < nsp; k++) { + resid[kspecial] -= CSoln[kins + k]; + } + } + } else { + kindexSP = 0; + for (isp = 0; isp < m_numSurfPhases; isp++) { + nsp = m_nSpeciesSurfPhase[isp]; + InterfaceKinetics *kinPtr = m_objects[isp]; + int surfIndex = kinPtr->surfacePhaseIndex(); + kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); + kins = kindexSP; + kinPtr->getNetProductionRates(DATA_PTR(m_netProductionRatesSave)); + for (k = 0; k < nsp; k++, kindexSP++) { + resid[kindexSP] = - m_netProductionRatesSave[kstart + k]; + } + kspecial = kins + m_spSurfLarge[isp]; + sd = m_ptrsSurfPhase[isp]->siteDensity(); + resid[kspecial] = sd; + for (k = 0; k < nsp; k++) { + resid[kspecial] -= CSoln[kins + k]; + } + } + } + + if (m_bulkFunc == BULK_DEPOSITION) { + kindexSP = m_numTotSurfSpecies; + for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + double *XBlk = DATA_PTR(m_numEqn1); + //ThermoPhase *THptr = m_bulkPhasePtrs[isp]; + //THptr->getMoleFractions(XBlk); + nsp = m_nSpeciesSurfPhase[isp]; + int surfPhaseIndex = m_indexKinObjSurfPhase[isp]; + InterfaceKinetics *m_kin = m_objects[isp]; + grRate = 0.0; + kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex); + for (k = 0; k < nsp; k++) { + if (m_netProductionRatesSave[kstart + k] > 0.0) { + grRate += m_netProductionRatesSave[kstart + k]; + } + } + resid[kindexSP] = m_bulkPhasePtrs[isp]->molarDensity(); + for (k = 0; k < nsp; k++) { + resid[kindexSP] -= CSoln[kindexSP + k]; + } + if (grRate > 0.0) { + for (k = 1; k < nsp; k++) { + if (m_netProductionRatesSave[kstart + k] > 0.0) { + resid[kindexSP + k] = XBlk[k] * grRate + - m_netProductionRatesSave[kstart + k]; + } else { + resid[kindexSP + k] = XBlk[k] * grRate; + } + } + } else { + grRate = 1.0E-6; + grRate += fabs(m_netProductionRatesSave[kstart + k]); + for (k = 1; k < nsp; k++) { + resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp); + } + } + if (do_time) { + for (k = 1; k < nsp; k++) { + resid[kindexSP + k] += + lenScale / deltaT * + (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]); + } + } + kindexSP += nsp; + } + } + } + } + + /* + * Calculate the Jacobian and residual + * + * @internal + * This routine uses the m_numEqn2 vector + * as temporary internal storage. + */ + void solveSP::resjac_eval(std::vector &JacCol, + double resid[], double CSoln[], + const double CSolnOld[], const bool do_time, + const double deltaT) + { + int kColIndex = 0, nsp, jsp, i, kCol; + double dc, cSave, sd; + double *col_j; + /* + * Calculate the residual + */ + fun_eval(resid, CSoln, CSolnOld, do_time, deltaT); + /* + * Now we will look over the columns perturbing each unknown. + */ + for (jsp = 0; jsp < m_numSurfPhases; jsp++) { + nsp = m_nSpeciesSurfPhase[jsp]; + sd = m_ptrsSurfPhase[jsp]->siteDensity(); + for (kCol = 0; kCol < nsp; kCol++) { + cSave = CSoln[kColIndex]; + dc = fmaxx(1.0E-10 * sd, fabs(cSave) * 1.0E-7); + CSoln[kColIndex] += dc; + fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT); + col_j = JacCol[kColIndex]; + for (i = 0; i < m_neq; i++) { + col_j[i] = (m_numEqn2[i] - resid[i])/dc; + } + CSoln[kColIndex] = cSave; + kColIndex++; + } + } + + if (m_bulkFunc == BULK_DEPOSITION) { + for (jsp = 0; jsp < m_numBulkPhasesSS; jsp++) { + nsp = m_numBulkSpecies[jsp]; + sd = m_bulkPhasePtrs[jsp]->molarDensity(); + for (kCol = 0; kCol < nsp; kCol++) { + cSave = CSoln[kColIndex]; + dc = fmaxx(1.0E-10 * sd, fabs(cSave) * 1.0E-7); + CSoln[kColIndex] += dc; + fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT); + col_j = JacCol[kColIndex]; + for (i = 0; i < m_neq; i++) { + col_j[i] = (m_numEqn2[i] - resid[i])/dc; + } + CSoln[kColIndex] = cSave; + kColIndex++; + } + } + } + } + + +#define APPROACH 0.80 + + static double calc_damping(double x[], double dxneg[], int dim, int *label) + + /* This function calculates a damping factor for the Newton iteration update + * vector, dxneg, to insure that all site and bulk fractions, x, remain + * bounded between zero and one. + * + * dxneg[] = negative of the update vector. + * + * The constant "APPROACH" sets the fraction of the distance to the boundary + * that the step can take. If the full step would not force any fraction + * outside of 0-1, then Newton's method is allowed to operate normally. + */ + + { + int i; + double damp = 1.0, xnew, xtop, xbot; + static double damp_old = 1.0; + + *label = -1; + + for (i = 0; i < dim; i++) { + + /* + * Calculate the new suggested new value of x[i] + */ + + xnew = x[i] - damp * dxneg[i]; + + /* + * Calculate the allowed maximum and minimum values of x[i] + * - Only going to allow x[i] to converge to zero by a + * single order of magnitude at a time + */ + + xtop = 1.0 - 0.1*fabs(1.0-x[i]); + xbot = fabs(x[i]*0.1) - 1.0e-16; + if (xnew > xtop ) { + damp = - APPROACH * (1.0 - x[i]) / dxneg[i]; + *label = i; + } + else if (xnew < xbot) { + damp = APPROACH * x[i] / dxneg[i]; + *label = i; + } else if (xnew > 3.0*MAX(x[i], 1.0E-10)) { + damp = - 2.0 * MAX(x[i], 1.0E-10) / dxneg[i]; + *label = i; + } + } + + if (damp < 1.0e-2) damp = 1.0e-2; + /* + * Only allow the damping parameter to increase by a factor of three each + * iteration. Heuristic to avoid oscillations in the value of damp + */ + + if (damp > damp_old*3) { + damp = damp_old*3; + *label = -1; + } + + /* + * Save old value of the damping parameter for use + * in subsequent calls. + */ + + damp_old = damp; + return damp; + + } /* calc_damping */ +#undef APPROACH + + /* + * This function calculates the norm of an update, dx[], + * based on the weighted values of x. + */ + static double calcWeightedNorm(const double wtX[], const double dx[], int dim) { + double norm = 0.0; + double tmp; + if (dim == 0) return 0.0; + for (int i = 0; i < dim; i++) { + tmp = dx[i] / wtX[i]; + norm += tmp * tmp; + } + return (sqrt(norm/dim)); + } + + /* + * Calculate the weighting factors for norms wrt both the species + * concentration unknowns and the residual unknowns. + * + */ + void solveSP::calcWeights(double wtSpecies[], double wtResid[], + const Array2D &Jac, const double CSoln[], + const double abstol, const double reltol) + { + int k, jcol, kindex, isp, nsp; + double sd; + /* + * First calculate the weighting factor for the concentrations of + * the surface species and bulk species. + */ + kindex = 0; + for (isp = 0; isp < m_numSurfPhases; isp++) { + nsp = m_nSpeciesSurfPhase[isp]; + sd = m_ptrsSurfPhase[isp]->siteDensity(); + for (k = 0; k < nsp; k++, kindex++) { + wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); + } + } + if (m_bulkFunc == BULK_DEPOSITION) { + for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + nsp = m_numBulkSpecies[isp]; + sd = m_bulkPhasePtrs[isp]->molarDensity(); + for (k = 0; k < nsp; k++, kindex++) { + wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); + } + } + } + /* + * Now do the residual Weights. Since we have the Jacobian, we + * will use it to generate a number based on the what a significant + * change in a solution variable does to each residual. + * This is a row sum scale operation. + */ + for (k = 0; k < m_neq; k++) { + wtResid[k] = 0.0; + for (jcol = 0; jcol < m_neq; jcol++) { + wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]); + } + } + } + + /* + * This routine calculates a pretty conservative 1/del_t based + * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees + * diagonal dominance. + * + * Small surface fractions are allowed to intervene in the del_t + * determination, no matter how small. This may be changed. + * Now minimum changed to 1.0e-12, + * + * Maximum time step set to time_scale. + */ + double solveSP:: + calc_t(double netProdRateSolnSP[], double XMolSolnSP[], + int *label, int *label_old, double *label_factor, int ioflag) + { + int k, isp, nsp, kstart; + double inv_timeScale = 1.0E-10; + double sden, tmp; + int kindexSP = 0; + *label = 0; + int ispSpecial = 0; + int kspSpecial = 0; + updateMFSolnSP(XMolSolnSP); + for (isp = 0; isp < m_numSurfPhases; isp++) { + nsp = m_nSpeciesSurfPhase[isp]; + + // Get the interface kinetics associated with this surface + InterfaceKinetics *m_kin = m_objects[isp]; + + // Calcuate the start of the species index for surfaces within + // the InterfaceKinetics object + int surfIndex = m_kin->surfacePhaseIndex(); + kstart = m_kin->kineticsSpeciesIndex(0, surfIndex); + ThermoPhase& THref = m_kin->thermo(surfIndex); + + m_kin->getNetProductionRates(DATA_PTR(m_numEqn1)); + + sden = THref.molarDensity(); + for (k = 0; k < nsp; k++, kindexSP++) { + int kspindex = kstart + k; + netProdRateSolnSP[kindexSP] = m_numEqn1[kspindex]; + if (XMolSolnSP[kindexSP] <= 1.0E-10) { + tmp = 1.0E-10; + } else { + tmp = XMolSolnSP[kindexSP]; + } + tmp *= sden; + tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp); + if (netProdRateSolnSP[kindexSP]> 0.0) tmp /= 100.; + if (tmp > inv_timeScale) { + inv_timeScale = tmp; + *label = kindexSP; + ispSpecial = isp; + kspSpecial = k; + } + } + } + + /* + * Increase time step exponentially as same species repeatedly + * controls time step + */ + if (*label == *label_old) { + *label_factor *= 1.5; + } else { + *label_old = *label; + *label_factor = 1.0; + } + inv_timeScale = inv_timeScale / *label_factor; +#ifdef DEBUG_SOLVESP + if (ioflag > 1) { + if (*label_factor > 1.0) { + printf("Delta_t increase due to repeated controlling species = %e\n", + *label_factor); + } + int kkin = m_kinSpecIndex[*label]; + InterfaceKinetics *m_kin = m_objects[ispSpecial]; + string sn = m_kin->kineticsSpeciesName(kkin); + printf("calc_t: spec=%d(%s) sf=%e pr=%e dt=%e\n", + *label, sn.c_str(), XMolSolnSP[*label], + netProdRateSolnSP[*label], 1.0/inv_timeScale); + } +#endif + + return (inv_timeScale); + + } /* calc_t */ + + + /** + * printResJac(): prints out the residual and Jacobian. + * + */ +#ifdef DEBUG_SOLVESP + void solveSP::printResJac(int ioflag, int neq, const Array2D &Jac, + double resid[], double wtRes[], + double norm) + { + int i, j, isp, nsp, irowKSI, irowISP; + int kstartKSI; + int kindexSP = 0; + string sname, pname, cname; + if (ioflag > 1) { + printf(" Printout of residual and jacobian\n"); + printf("\t Residual: weighted norm = %10.4e\n", norm); + printf("\t Index Species_Name Residual " + "Resid/wtRes wtRes\n"); + for (isp = 0; isp < m_numSurfPhases; isp++) { + nsp = m_nSpeciesSurfPhase[isp]; + InterfaceKinetics *m_kin = m_objects[isp]; + int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; + m_kin->getNetProductionRates(DATA_PTR(m_numEqn1)); + kstartKSI = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex); + SurfPhase *sp_ptr = m_ptrsSurfPhase[isp]; + pname = sp_ptr->id(); + + for (int k = 0; k < nsp; k++, kindexSP++) { + sname = sp_ptr->speciesName(k); + cname = pname + ":" + sname; + printf("\t %d: %-24s: %11.3e %11.3e %11.3e\n", kindexSP, + cname.c_str(), resid[kindexSP], + resid[kindexSP]/wtRes[kindexSP], wtRes[kindexSP]); + } + } + if (m_bulkFunc == BULK_DEPOSITION) { + for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + // fill in + } + } + if (ioflag > 2) { + printf("\t Jacobian:\n"); + for (i = 0; i < m_neq; i++) { + irowISP = m_kinObjIndex[i]; + InterfaceKinetics *m_kin = m_objects[irowISP]; + irowKSI = m_kinSpecIndex[i]; + ThermoPhase& THref = m_kin->speciesPhase(irowKSI); + int phaseIndex = m_kin->speciesPhaseIndex(irowKSI); + kstartKSI = m_kin->kineticsSpeciesIndex(0, phaseIndex); + int klocal = i - m_eqnIndexStartSolnPhase[irowISP]; + sname = THref.speciesName(klocal); + printf("\t Row %d:%-16s:\n", i, sname.c_str()); + printf("\t "); + for (j = 0; j < m_neq; j++) { + printf("%10.4e ", Jac(i,j)); + } + printf("\n"); + } + } + } + } /* printResJac */ +#endif + + /* + * Optional printing at the start of the solveSP problem + */ + void solveSP::print_header(int ioflag, int ifunc, double time_scale, + int damping, double reltol, double abstol, + double TKelvin, + double PGas, double netProdRate[], + double XMolKinSpecies[]) { + if (ioflag) { + printf("\n================================ SOLVESP CALL SETUP " + "========================================\n"); + if (ifunc == SFLUX_INITIALIZE) { + printf("\n SOLVESP Called with Initialization turned on\n"); + printf(" Time scale input = %9.3e\n", time_scale); + } + else if (ifunc == SFLUX_RESIDUAL) { + printf("\n SOLVESP Called to calculate steady state residual\n"); + printf( " from a good initial guess\n"); + } + else if (ifunc == SFLUX_JACOBIAN) { + printf("\n SOLVESP Called to calculate steady state jacobian\n"); + printf( " from a good initial guess\n"); + } + else if (ifunc == SFLUX_TRANSIENT) { + printf("\n SOLVESP Called to integrate surface in time\n"); + printf( " for a total of %9.3e sec\n", time_scale); + } + else { + fprintf(stderr,"Unknown ifunc flag = %d\n", ifunc); + exit (-1); + } + + if (m_bulkFunc == BULK_DEPOSITION) + printf(" The composition of the Bulk Phases will be calculated\n"); + else if (m_bulkFunc == BULK_ETCH) + printf(" Bulk Phases have fixed compositions\n"); + else { + fprintf(stderr,"Unknown bulkFunc flag = %d\n", m_bulkFunc); + exit (-1); + } + + if (damping) + printf(" Damping is ON \n"); + else + printf(" Damping is OFF \n"); + + printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol); + } + + /* + * Print out the initial guess + */ +#ifdef DEBUG_SOLVESP + if (ioflag > 1) { + printf("\n================================ INITIAL GUESS " + "========================================\n"); + int kindexSP = 0; + for (int isp = 0; isp < m_numSurfPhases; isp++) { + InterfaceKinetics *m_kin = m_objects[isp]; + int surfIndex = m_kin->surfacePhaseIndex(); + int nPhases = m_kin->nPhases(); + m_kin->getNetProductionRates(netProdRate); + updateMFKinSpecies(XMolKinSpecies, isp); + + printf("\n IntefaceKinetics Object # %d\n\n", isp); + + printf("\t Number of Phases = %d\n", nPhases); + printf("\t Temperature = %10.3e Kelvin\n", TKelvin); + printf("\t Pressure = %10.3g Pa\n\n", PGas); + printf("\t Phase:SpecName Prod_Rate MoleFraction kindexSP\n"); + printf("\t -------------------------------------------------------" + "----------\n"); + + int kspindex = 0; + bool inSurfacePhase = false; + for (int ip = 0; ip < nPhases; ip++) { + if (ip == surfIndex) { + inSurfacePhase = true; + } else { + inSurfacePhase = false; + } + ThermoPhase &THref = m_kin->thermo(ip); + int nsp = THref.nSpecies(); + string pname = THref.id(); + for (int k = 0; k < nsp; k++) { + string sname = THref.speciesName(k); + string cname = pname + ":" + sname; + if (inSurfacePhase) { + printf("\t %-24s %10.3e %10.3e %d\n", cname.c_str(), + netProdRate[kspindex], XMolKinSpecies[kspindex], + kindexSP); + kindexSP++; + } else { + printf("\t %-24s %10.3e %10.3e\n", cname.c_str(), + netProdRate[kspindex], XMolKinSpecies[kspindex]); + } + kspindex++; + } + } + printf("==========================================================" + "=================================\n"); + } + } +#endif + if (ioflag == 1) { + printf("\n\n\t Iter Time Del_t Damp DelX " + " Resid Name-Time Name-Damp\n"); + printf( "\t -----------------------------------------------" + "------------------------------------\n"); + } + } + + void solveSP::printIteration(int ioflag, double damp, int label_d, + int label_t, + double inv_t, double t_real, int iter, + double update_norm, double resid_norm, + double netProdRate[], double CSolnSP[], + double resid[], double XMolSolnSP[], + double wtSpecies[], int dim, bool do_time) + { + int i, k; + string nm; + if (ioflag == 1) { + + printf("\t%6d ", iter); + if (do_time) + printf("%9.4e %9.4e ", t_real, 1.0/inv_t); + else + for (i = 0; i < 22; i++) printf(" "); + if (damp < 1.0) + printf("%9.4e ", damp); + else + for (i = 0; i < 11; i++) printf(" "); + printf("%9.4e %9.4e", update_norm, resid_norm); + if (do_time) { + k = m_kinSpecIndex[label_t]; + int isp = m_kinObjIndex[label_t]; + InterfaceKinetics *m_kin = m_objects[isp]; + nm = m_kin->kineticsSpeciesName(k); + printf(" %-16s", nm.c_str()); + } else { + for (i = 0; i < 16; i++) printf(" "); + } + if (label_d >= 0) { + k = m_kinSpecIndex[label_d]; + int isp = m_kinObjIndex[label_d]; + InterfaceKinetics *m_kin = m_objects[isp]; + nm = m_kin->kineticsSpeciesName(k); + printf(" %-16s", nm.c_str()); + } + printf("\n"); + } +#ifdef DEBUG_SOLVESP + else if (ioflag > 1) { + + updateMFSolnSP(XMolSolnSP); + printf("\n\t Weighted norm of update = %10.4e\n", update_norm); + + printf("\t Name Prod_Rate XMol Conc " + " Conc_Old wtConc"); + if (damp < 1.0) printf(" UnDamped_Conc"); + printf("\n"); + printf("\t---------------------------------------------------------" + "-----------------------------\n"); + int kindexSP = 0; + for (int isp = 0; isp < m_numSurfPhases; isp++) { + int nsp = m_nSpeciesSurfPhase[isp]; + InterfaceKinetics *m_kin = m_objects[isp]; + //int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; + m_kin->getNetProductionRates(DATA_PTR(m_numEqn1)); + for (int k = 0; k < nsp; k++, kindexSP++) { + int kspIndex = m_kinSpecIndex[kindexSP]; + nm = m_kin->kineticsSpeciesName(kspIndex); + printf("\t%-16s %10.3e %10.3e %10.3e %10.3e %10.3e ", + nm.c_str(), + m_numEqn1[kspIndex], + XMolSolnSP[kindexSP], + CSolnSP[kindexSP], CSolnSP[kindexSP]+damp*resid[kindexSP], + wtSpecies[kindexSP]); + if (damp < 1.0) { + printf("%10.4e ", CSolnSP[kindexSP]+(damp-1.0)*resid[kindexSP]); + if (label_d == kindexSP) printf(" Damp "); + } + if (label_t == kindexSP) printf(" Tctrl"); + printf("\n"); + } + + } + + printf("\t--------------------------------------------------------" + "------------------------------\n"); + } +#endif + } /* printIteration */ + + + void solveSP::printFinal(int ioflag, double damp, int label_d, int label_t, + double inv_t, double t_real, int iter, + double update_norm, double resid_norm, + double netProdRateKinSpecies[], const double CSolnSP[], + const double resid[], double XMolSolnSP[], + const double wtSpecies[], const double wtRes[], + int dim, bool do_time, + double TKelvin, double PGas) + { + int i, k; + string nm; + if (ioflag == 1) { + + printf("\tFIN%3d ", iter); + if (do_time) + printf("%9.4e %9.4e ", t_real, 1.0/inv_t); + else + for (i = 0; i < 22; i++) printf(" "); + if (damp < 1.0) + printf("%9.4e ", damp); + else + for (i = 0; i < 11; i++) printf(" "); + printf("%9.4e %9.4e", update_norm, resid_norm); + if (do_time) { + k = m_kinSpecIndex[label_t]; + int isp = m_kinObjIndex[label_t]; + InterfaceKinetics *m_kin = m_objects[isp]; + nm = m_kin->kineticsSpeciesName(k); + printf(" %-16s", nm.c_str()); + } else { + for (i = 0; i < 16; i++) printf(" "); + } + if (label_d >= 0) { + k = m_kinSpecIndex[label_d]; + int isp = m_kinObjIndex[label_d]; + InterfaceKinetics *m_kin = m_objects[isp]; + nm = m_kin->kineticsSpeciesName(k); + printf(" %-16s", nm.c_str()); + } + printf(" -- success\n"); + } +#ifdef DEBUG_SOLVESP + else if (ioflag > 1) { + + + printf("\n================================== FINAL RESULT =========" + "==================================================\n"); + updateMFSolnSP(XMolSolnSP); + printf("\n Weighted norm of solution update = %10.4e\n", update_norm); + printf(" Weighted norm of residual update = %10.4e\n\n", resid_norm); + + printf(" Name Prod_Rate XMol Conc " + " wtConc Resid Resid/wtResid wtResid"); + if (damp < 1.0) printf(" UnDamped_Conc"); + printf("\n"); + printf("---------------------------------------------------------------" + "---------------------------------------------\n"); + int kindexSP = 0; + for (int isp = 0; isp < m_numSurfPhases; isp++) { + int nsp = m_nSpeciesSurfPhase[isp]; + InterfaceKinetics *m_kin = m_objects[isp]; + //int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; + m_kin->getNetProductionRates(DATA_PTR(m_numEqn1)); + for (int k = 0; k < nsp; k++, kindexSP++) { + int kspIndex = m_kinSpecIndex[kindexSP]; + nm = m_kin->kineticsSpeciesName(kspIndex); + printf("%-16s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e", + nm.c_str(), + m_numEqn1[kspIndex], + XMolSolnSP[kindexSP], + CSolnSP[kindexSP], + wtSpecies[kindexSP], + resid[kindexSP], + resid[kindexSP]/wtRes[kindexSP], wtRes[kindexSP]); + if (damp < 1.0) { + printf("%10.4e ", CSolnSP[kindexSP]+(damp-1.0)*resid[kindexSP]); + if (label_d == kindexSP) printf(" Damp "); + } + if (label_t == kindexSP) printf(" Tctrl"); + printf("\n"); + } + + } + printf("---------------------------------------------------------------" + "---------------------------------------------\n"); + double *XMolKinSpecies = DATA_PTR(m_numEqn2); + kindexSP = 0; + for (int isp = 0; isp < m_numSurfPhases; isp++) { + InterfaceKinetics *m_kin = m_objects[isp]; + int surfIndex = m_kin->surfacePhaseIndex(); + int nPhases = m_kin->nPhases(); + m_kin->getNetProductionRates(netProdRateKinSpecies); + + updateMFKinSpecies(XMolKinSpecies, isp); + + printf("\n IntefaceKinetics Object # %d\n\n", isp); + + printf("\t Number of Phases = %d\n", nPhases); + printf("\t Temperature = %10.3e Kelvin\n", TKelvin); + printf("\t Pressure = %10.3g Pa\n\n", PGas); + printf("\t Phase:SpecName Prod_Rate MoleFraction kindexSP\n"); + printf("\t--------------------------------------------------------------" + "---\n"); + + int kspindex = 0; + bool inSurfacePhase = false; + for (int ip = 0; ip < nPhases; ip++) { + if (ip == surfIndex) { + inSurfacePhase = true; + } else { + inSurfacePhase = false; + } + ThermoPhase &THref = m_kin->thermo(ip); + int nsp = THref.nSpecies(); + string pname = THref.id(); + for (k = 0; k < nsp; k++) { + string sname = THref.speciesName(k); + string cname = pname + ":" + sname; + if (inSurfacePhase) { + printf("\t%-24s %10.3e %10.3e %d\n", cname.c_str(), + netProdRateKinSpecies[kspindex], XMolKinSpecies[kspindex], kindexSP); + kindexSP++; + } else { + printf("\t%-24s %10.3e %10.3e\n", cname.c_str(), + netProdRateKinSpecies[kspindex], XMolKinSpecies[kspindex]); + } + kspindex++; + } + } + } + printf("\n"); + printf("===============================================================" + "============================================\n\n"); + } +#endif + } + +#ifdef DEBUG_SOLVESP + void solveSP:: + printIterationHeader(int ioflag, double damp,double inv_t, double t_real, + int iter, bool do_time) + { + if (ioflag > 1) { + printf("\n===============================Iteration %5d " + "=================================\n", iter); + if (do_time) { + printf(" Transient step with: Real Time_n-1 = %10.4e sec,", t_real); + printf(" Time_n = %10.4e sec\n", t_real + 1.0/inv_t); + printf(" Delta t = %10.4e sec", 1.0/inv_t); + } else { + printf(" Steady Solve "); + } + if (damp < 1.0) { + printf(", Damping value = %10.4e\n", damp); + } else { + printf("\n"); + } + } + } +#endif + +} diff --git a/Cantera/src/kinetics/solveSP.h b/Cantera/src/kinetics/solveSP.h new file mode 100644 index 000000000..c4578cde5 --- /dev/null +++ b/Cantera/src/kinetics/solveSP.h @@ -0,0 +1,686 @@ +/** + * @file solveSP.h + * Header file for implicit surface problem solver + * (see \ref kinetics and class \link Cantera::solveSP solveSP\endlink). + */ +/* + * $Id$ + */ +/* + * Copywrite 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + * See file License.txt for licensing information. + */ + +#ifndef SOLVESP_H +#define SOLVESP_H + +#include "ImplicitSurfChem.h" +#include "InterfaceKinetics.h" + +#include +#include "Array.h" + +//! Solution Methods +/*! + * Flag to specify the solution method + * + * 1: SFLUX_INITIALIZE = This assumes that the initial guess supplied to the + * routine is far from the correct one. Substantial + * work plus transient time-stepping is to be expected + * to find a solution. + * 2: SFLUX_RESIDUAL = Need to solve the surface problem in order to + * calculate the surface fluxes of gas-phase species. + * (Can expect a moderate change in the solution + * vector -> try to solve the system by direct + * methods + * with no damping first -> then, try time-stepping + * if the first method fails) + * A "time_scale" supplied here is used in the + * algorithm to determine when to shut off + * time-stepping. + * 3: SFLUX_JACOBIAN = Calculation of the surface problem is due to the + * need for a numerical jacobian for the gas-problem. + * The solution is expected to be very close to the + * initial guess, and accuracy is needed. + * 4: SFLUX_TRANSIENT = The transient calculation is performed here for an + * amount of time specified by "time_scale". It is + * not garraunted to be time-accurate - just stable + * and fairly fast. The solution after del_t time is + * returned, whether it's converged to a steady + * state or not. + */ +const int SFLUX_INITIALIZE = 1; +const int SFLUX_RESIDUAL = 2; +const int SFLUX_JACOBIAN = 3; +const int SFLUX_TRANSIENT = 4; + + +/* + * bulkFunc: Functionality expected from the bulk phase. This changes the + * equations that will be used to solve for the bulk mole + * fractions. + * 1: BULK_DEPOSITION = deposition of a bulk phase is to be expected. + * Bulk mole fractions are determined from ratios of + * growth rates of bulk species. + * 2: BULK_ETCH = Etching of a bulk phase is to be expected. + * Bulk mole fractions are assumed constant, and given + * by the initial conditions. This is also used + whenever the condensed phase is part of the larger + solution. +*/ +const int BULK_DEPOSITION = 1; +const int BULK_ETCH = 2; + +namespace Cantera { + + class InterfaceKinetics; + + //! Method to solve a pseudo steady state surface problem + /*! + * The following class handles solving the surface problem. + * The calculation uses Newton's method to + * obtain the surface fractions of the surface and bulk species by + * requiring that the + * surface species production rate = 0 and that the either the + * bulk fractions are proportional to their production rates + * or they are constants. + * + * Currently, the bulk mole fractions are treated as constants. + * Implementation of their being added to the unknown solution + * vector is delayed. + * + * Lets introduce the unknown vector for the "surface + * problem". The surface problem is defined as the evaluation of the surface + * site fractions for multiple surface phases. + * The unknown vector will consist of the vector of surface concentrations for each + * species in each surface vector. Species are grouped first by their surface phases + * + * C_i_j = Concentration of the ith species in the jth surface phase + * Nj = number of surface species in the jth surface phase + * + * The unknown solution vector is defined as follows: + * + * kindexSP + * ---------------------------- + * C_0_0 0 + * C_1_0 1 + * C_2_0 2 + * . . . ... + * C_N0-1_0 N0-1 + * C_0_1 N0 + * C_1_1 N0+1 + * C_2_1 N0+2 + * . . . ... + * C_N1-1_1 NO+N1-1 + * + * + * Note there are a couple of different types of species indecices + * floating around in the formulation of this object. + * + * kindexSP This is the species index in the contiguous vector of unknowns + * for the surface problem. + * + * Note, in the future, BULK_DEPOSITION systems will be added, and the solveSP unknown + * vector will get more complicated. It will include the mole fraction and growth rates + * of specified bulk phases + * + * Indecises which relate to individual kinetics objects use the suffix KSI (kinetics + * species index). + * + * + * Solution Method + * + * This routine is typically used within a residual calculation in a large code. + * It's typically invoked millions of times for large calculations, and it must + * work every time. Therefore, requirements demand that it be robust but also + * efficient. + * + * The solution methodology is largely determined by the ifunc<\TT> parameter, + * that is input to the solution object. This parameter may have the following + * 4 values: + * + * + * 1: SFLUX_INITIALIZE = This assumes that the initial guess supplied to the + * routine is far from the correct one. Substantial + * work plus transient time-stepping is to be expected + * to find a solution. + * + * 2: SFLUX_RESIDUAL = Need to solve the surface problem in order to + * calculate the surface fluxes of gas-phase species. + * (Can expect a moderate change in the solution + * vector -> try to solve the system by direct methods + * with no damping first -> then, try time-stepping + * if the first method fails) + * A "time_scale" supplied here is used in the + * algorithm to determine when to shut off + * time-stepping. + * + * 3: SFLUX_JACOBIAN = Calculation of the surface problem is due to the + * need for a numerical jacobian for the gas-problem. + * The solution is expected to be very close to the + * initial guess, and extra accuracy is needed because + * solution variables have been delta'd from + * nominal values to create jacobian entries. + * + * 4: SFLUX_TRANSIENT = The transient calculation is performed here for an + * amount of time specified by "time_scale". It is + * not garraunted to be time-accurate - just stable + * and fairly fast. The solution after del_t time is + * returned, whether it's converged to a steady + * state or not. This is a poor man's time stepping + * algorithm. + * + * Psuedo time stepping algorithm: + * The time step is determined from sdot[], so so that the time step + * doesn't ever change the value of a variable by more than 100%. + * + * This algorithm does use a damped Newton's method to relax the equations. + * Damping is based on a "delta damping" technique. The solution unknowns + * are not allowed to vary too much between iterations. + * + * + * EXTRA_ACCURACY:A constant that is the ratio of the required update norm in + * this Newton iteration compared to that in the nonlinear solver. + * A value of 0.1 is used so surface species are safely overconverged. + * + * Functions called: + *---------------------------------------------------------------------------- + * + * ct_dgetrf -- First half of LAPACK direct solve of a full Matrix + * + * ct_dgetrs -- Second half of LAPACK direct solve of a full matrix. Returns + * solution vector in the right-hand-side vector, resid. + * + *---------------------------------------------------------------------------- + * + */ + class solveSP { + + public: + + //! Constructor for the object + /*! + * @param surfChemPtr Pointer to the ImplicitSurfChem object that + * defines the surface problem to be solved. + * + * @param bulkFunc Integer representing how the bulk phases + * should be handled. Currently, only the + * default value of BULK_ETCH is supported. + */ + solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc = BULK_ETCH); + + //! Destructor. Deletes the integrator. + ~solveSP(); + + private: + + //! Unimplemented private copy constructor + solveSP(const solveSP &right); + + //! Unimplemented private assignment operator + solveSP& operator=(const solveSP &right); + + public: + + //! Main routine that actually calculates the pseudo steady state + //! of the surface problem + /*! + * The actual converged solution is returned as part of the + * internal state of the InterfaceKinetics objects. + * + * @param ifunc Determines the type of solution algorithm to be + * used. Possible values are SFLUX_INITIALIZE , + * SFLUX_RESIDUAL SFLUX_JACOBIAN SFLUX_TRANSIENT . + * + * @param time_scale Time over which to integrate the surface equations, + * where applicable + * + * @param TKelvin Temperature (kelvin) + * + * @param PGas Pressure (pascals) + * + * @param reltol Relative tolerance to use + * @param abstol absolute tolerance. + * + * @return Returns 1 if the surface problem is successfully solved. + * Returns -1 if the surface problem wasn't solved successfully. + * Note the actual converged solution is returned as part of the + * internal state of the InterfaceKinetics objects. + */ + int solveSurfProb(int ifunc, double time_scale, double TKelvin, + double PGas, double reltol, double abstol); + + private: + + //! Printing routine that gets called at the start of every + //! invocation + void print_header(int ioflag, int ifunc, double time_scale, + int damping, double reltol, double abstol, + double TKelvin, double PGas, double netProdRate[], + double XMolKinSpecies[]); + +#ifdef DEBUG_SOLVESP + + void printResJac(int ioflag, int neq, const Array2D &Jac, + double resid[], double wtResid[], double norm); +#endif + + //! Printing routine that gets called after every iteration + void printIteration(int ioflag, double damp, int label_d, int label_t, + double inv_t, double t_real, int iter, + double update_norm, double resid_norm, + double netProdRate[], double CSolnSP[], + double resid[], double XMolSolnSP[], + double wtSpecies[], int dim, bool do_time); + + + //! Print a summary of the solution + /*! + * + */ + void printFinal(int ioflag, double damp, int label_d, int label_t, + double inv_t, double t_real, int iter, + double update_norm, double resid_norm, + double netProdRateKinSpecies[], const double CSolnSP[], + const double resid[], double XMolSolnSP[], + const double wtSpecies[], const double wtRes[], + int dim, bool do_time, + double TKelvin, double PGas); + + //! Calculate a conservative delta T to use in a pseudo-steady state + //! algorithm + /*! + * This routine calculates a pretty conservative 1/del_t based + * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees + * diagonal dominance. + * + * Small surface fractions are allowed to intervene in the del_t + * determination, no matter how small. This may be changed. + * Now minimum changed to 1.0e-12, + * + * Maximum time step set to time_scale. + * + * @param netProdRateSolnSP Output variable. Net production rate + * of all of the species in the solution vector. + * @param XMolSolnSP output variable. + * Mole fraction of all of the species in the solution vector + * @param label Output variable. Pointer to the value of the + * species index (kindexSP) that is controlling + * the time step + * @param label_old Output variable. Pointer to the value of the + * species index (kindexSP) that controlled + * the time step at the previous iteration + * @param label_factor Output variable. Pointer to the current + * factor that is used to indicate the same species + * is controlling the time step. + * + * @param ioflag Level of the output requested. + * + * @return Returns the 1. / delta T to be used on the next step + */ + double calc_t(double netProdRateSolnSP[], double XMolSolnSP[], + int *label, int *label_old, + double *label_factor, int ioflag); + + //! Calculate the solution and residual weights + /*! + * @param wtSpecies Weights to use for the soln unknowns. These + * are in concentration units + * @param wtResid Weights to sue for the residual unknowns. + * + * @param Jac Jacobian. Row sum scaling is used for the Jacobian + * @param CSolnSP Solution vector for the surface problem + * @param abstol Absolute error tolerance + * @param reltol Relative error tolerance + */ + void calcWeights(double wtSpecies[], double wtResid[], + const Array2D &Jac, const double CSolnSP[], + const double abstol, const double reltol); + +#ifdef DEBUG_SOLVESP + //! Utility routine to print a header for high lvls of debugging + /*! + * @param ioflag Lvl of debugging + * @param damp lvl of damping + * @param inv_t Inverse of the value of delta T + * @param t_real Value of the time + * @param iter Interation number + * @param do_time boolean indicating whether time stepping is taking + * place + */ + void printIterationHeader(int ioflag, double damp, + double inv_t, double t_real, int iter, + bool do_time); +#endif + + /** + * Update the surface states of the surface phases. + */ + void updateState(const double* cSurfSpec); + + //! Update mole fraction vector consisting of unknowns in surface problem + /*! + * @param XMolSolnSP Vector of mole fractions for the unknowns in the + * surface problem. + */ + void updateMFSolnSP(double *XMolSolnSP); + + //! Update the mole fraction vector for a specific kinetic species vector + //! corresponding to one InterfaceKinetics object + /*! + * @param XMolKinSp Mole fraction vector corresponding to a particular + * kinetic species for a single InterfaceKinetics Object + * This is a vector over all the species in all of the + * phases in the InterfaceKinetics object + * @param isp ID of the InterfaceKinetics Object. + */ + void updateMFKinSpecies(double *XMolKinSp, int isp); + + + //! Update the vector that keeps track of the largest species in each + //! surface phase. + /*! + * @param CsolnSP Vector of the current values of the surface concentrations + * in all of the surface species. + */ + void evalSurfLarge(const double *CSolnSP); + + //! Main Function evalulation + /*! + * + * @param resid output Vector of residuals, length = m_neq + * @param CSolnSP Vector of species concentrations, unknowns in the + * problem, length = m_neq + * @param CSolnSPOld Old Vector of species concentrations, unknowns in the + * problem, length = m_neq + * @param do_time Calculate a time dependent residual + * @param deltaT Delta time for time dependent problem. + */ + void fun_eval(double* resid, const double *CSolnSP, + const double *CSolnOldSP, const bool do_time, const double deltaT); + + //! Main routine that calculates the current residual and Jacobian + /*! + * @param JacCol Vector of pointers to the tops of columns of the + * Jacobian to be evalulated. + * @param resid output Vector of residuals, length = m_neq + * @param CSolnSP Vector of species concentrations, unknowns in the + * problem, length = m_neq. These are tweaked in order + * to derive the columns of the jacobian. + * @param CSolnSPOld Old Vector of species concentrations, unknowns in the + * problem, length = m_neq + * @param do_time Calculate a time dependent residual + * @param deltaT Delta time for time dependent problem. + */ + void resjac_eval(std::vector& JacCol, double* resid, + double *CSolnSP, + const double *CSolnSPOld, const bool do_time, + const double deltaT); + + //! Pointer to the manager of the implicit surface chemistry + //! problem + /*! + * This object actually calls the current object. Thus, we are + * providing a loop-back functionality here. + */ + ImplicitSurfChem *m_SurfChemPtr; + + //! Vector of interface kinetics objects + /*! + * Each of these is associated with one and only one surface phase. + */ + std::vector &m_objects; + + //! Total number of equations to solve in the implicit problem. + /*! + * Note, this can be zero, and frequently is + */ + int m_neq; + + //! This variable determines how the bulk phases are to be handled + /*! + * = BULK_ETCH (default) The concentrations of the bulk phases are + * considered constant, just as the gas phase is. + * They are not part of the solution vector. + * = BULK_DEPOSITION = + * We solve here for the composition of the bulk + * phases by calculating a growth rate. The equations + * for the species in the bulk phases are + * unknowns in this calculation. + */ + int m_bulkFunc; + + //! Number of surface phases in the surface problem + /*! + * This number is equal to the number of InterfaceKinetics objects + * in the problem. (until further noted) + */ + int m_numSurfPhases; + + //! Total number of surface species in all surface phases. + /*! + * This is also the number of equations to solve for m_mode=0 system + * It's equal to the sum of the number of species in each of the + * m_numSurfPhases. + */ + int m_numTotSurfSpecies; + + //! Mapping between the surface phases and the InterfaceKinetics objects + /*! + * Currently this is defined to be a 1-1 mapping (and probably assumed + * in some places) + * m_surfKinObjID[i] = i + */ + std::vector m_indexKinObjSurfPhase; + + //! Vector of length number of surface phases containing + //! the number of surface species in each phase + /*! + * Length is equal to the number of surface phases, m_numSurfPhases + */ + std::vector m_nSpeciesSurfPhase; + + //! Vector of surface phase pointers + /*! + * This is created during the constructor + * Length is equal to the number of surface phases, m_numSurfPhases + */ + std::vector m_ptrsSurfPhase; + + //! Index of the start of the unknowns for each solution phase + /*! + * i_eqn = m_eqnIndexStartPhase[isp] + * + * isp is the phase id in the list of phases solved by the + * surface problem. + * + * i_eqn is the equation number of the first unknown in the + * solution vector corresponding to isp'th phase. + */ + std::vector m_eqnIndexStartSolnPhase; + + //! Phase ID in the InterfaceKinetics object of the surface phase + /*! + * For each surface phase, this lists the PhaseId of the + * surface phase in the corresponding InterfaceKinetics object + * + * Length is equal to m_numSurfPhases + */ + std::vector m_kinObjPhaseIDSurfPhase; + + //! Total number of volumetric condensed phases included in the steady state + //! problem handled by this routine. + /*! + * This is equal to or less + * than the total number of volumetric phases in all of the InterfaceKinetics + * objects. We usually do not include bulk phases. Bulk phases + * are only included in the calculation when their domain isn't included + * in the underlying continuum model conservation equation system. + * + * This is equal to 0, for the time being + */ + int m_numBulkPhasesSS; + + //! Vector of number of species in the m_numBulkPhases phases. + /*! + * Length is number of bulk phases + */ + std::vector m_numBulkSpecies; + + //std::vector m_bulkKinObjID; + //std::vector m_bulkKinObjPhaseID; + + + //! Total number of species in all bulk phases. + /*! + * This is also the number of bulk equations to solve when bulk + * equation solving is turned on. + */ + int m_numTotBulkSpeciesSS; + + //! Vector of bulk phase pointers, length is equal to m_numBulkPhases. + /*! + * + */ + std::vector m_bulkPhasePtrs; + + //! Index between the equation index and the position in the + //! kinetic species array for the appropriate kinetics + //! operator + /*! + * Length = m_neq. + * + * ksp = m_kinSpecIndex[ieq] + * ksp is the kinetic species index for the ieq'th equation. + */ + std::vector m_kinSpecIndex; + + //! Index between the equation index and the index of the + //! InterfaceKinetics object + /*! + * Length m_neq + */ + std::vector m_kinObjIndex; + + //! Vector containing the indecies of the largest species + //! in each surface phase + /*! + * k = m_spSurfLarge[i] + * where + * k is the local species index, i.e., + * it varies from 0 num species in phase-1 + * i is the surface phase index in the problem + * + * length is equal to m_numSurfPhases + */ + std::vector m_spSurfLarge; + + //! m_atol is the absolute tolerance in real units. + /*! + * units are (kmol/m2) + */ + double m_atol; + + //! m_rtol is the relative error tolerance. + double m_rtol; + + //! maximum value of the time step + /*! + * units = seconds + */ + doublereal m_maxstep; + + //! Maximum number of species in any single kinetics operator + //! -> also maxed wrt the total # of solution species + int m_maxTotSpecies; + + //! Temporary vector with length equal to max m_maxTotSpecies + vector_fp m_netProductionRatesSave; + + //! Temporary vector with length equal to max m_maxTotSpecies + vector_fp m_numEqn1; + + //! Temporary vector with length equal to max m_maxTotSpecies + vector_fp m_numEqn2; + + //! Temporary vector with length equal to max m_maxTotSpecies + vector_fp m_CSolnSave; + + //! Solution vector + /*! + * length MAX(1, m_neq) + */ + vector_fp m_CSolnSP; + + //! Saved inital solution vector + /*! + * length MAX(1, m_neq) + */ + vector_fp m_CSolnSPInit; + + //! Saved solution vector at the old time step + /*! + * length MAX(1, m_neq) + */ + vector_fp m_CSolnSPOld; + + //! Weights for the residual norm calculation + /*! + * length MAX(1, m_neq) + */ + vector_fp m_wtResid; + + //! Weights for the species concentrations norm calculation + /*! + * length MAX(1, m_neq) + */ + vector_fp m_wtSpecies; + + //! Residual for the surface problem + /*! + * The residual vector of length "dim" that, that has the value + * of "sdot" for surface species. The residuals for the bulk + * species are a function of the sdots for all species in the bulk + * phase. The last residual of each phase enforces {Sum(fractions) + * = 1}. After linear solve (dgetrf_ & dgetrs_), resid holds the + * update vector. + * + * length MAX(1, m_neq) + */ + vector_fp m_resid; + + //! Vector of mole fractions + /*! + *length m_maxTotSpecies + */ + vector_fp m_XMolKinSpecies; + + //! pivots + /*! + * length MAX(1, m_neq) + */ + vector_int m_ipiv; + + //! Vector of pointers to the top of the columns of the + //! jacobians + /*! + * The "dim" by "dim" computed Jacobian matrix for the + * local Newton's method. + */ + std::vector m_JacCol; + + //! Jacobian + /*! + * m_neq by m_neq computed Jacobian matrix for the + * local Newton's method. + */ + Array2D m_Jac; + + + public: + int ioflag; + }; +} +#endif