GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2026-08-02 01:27:57 -05:00
Code Issues Packages Projects Releases Wiki Activity

Remove unused ImplicitSurfChem and SolveSP classes

Browse Source
This commit is contained in:
Ray Speth
2026-01-31 08:41:16 -06:00
committed by Ingmar Schoegl
parent 0abf617b6f
commit bf4554d4e3
9 changed files with 1 additions and 2113 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/cantera/kinetics/ImplicitSurfChem.h
-304
View File
@@ -1,304 +0,0 @@
/**
* @file ImplicitSurfChem.h
* Declarations for the implicit integration of surface site density equations
* (see @ref kineticsmgr and class
* @link Cantera::ImplicitSurfChem ImplicitSurfChem@endlink).
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef CT_IMPSURFCHEM_H
#define CT_IMPSURFCHEM_H
#include "cantera/numerics/Integrator.h"
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/kinetics/solveSP.h"
namespace Cantera
{
//! @defgroup surfSolverGroup Surface Problem Solver
//! Advances the surface coverages of the associated set of SurfacePhase
//! objects in time
/*!
* This function advances a set of SurfPhase objects, each associated with one
* InterfaceKinetics object, in time. The following equation is used for each
* surface phase, *i*.
*
* @f[
* \dot \theta_k = \dot s_k (\sigma_k / s_0)
* @f]
*
* In this equation,
* - @f$ \theta_k @f$ is the site coverage for the kth species.
* - @f$ \dot s_k @f$ is the source term for the kth species
* - @f$ \sigma_k @f$ is the number of surface sites covered by each species k.
* - @f$ s_0 @f$ is the total site density of the interfacial phase.
*
* Additionally, the 0'th equation in the set is discarded. Instead the
* alternate equation is solved for
*
* @f[
* \sum_{k=0}^{N-1} \dot \theta_k = 0
* @f]
*
* This last equation serves to ensure that sum of the @f$ \theta_k @f$ values
* stays constant.
*
* The object uses the CVODE software to advance the surface equations.
*
* The solution vector used by this object is as follows: For each surface
* phase with @f$ N_s @f$ surface sites, it consists of the surface coverages
* @f$ \theta_k @f$ for @f$ k = 0, N_s - 1 @f$
*
* @ingroup surfSolverGroup
*/
class ImplicitSurfChem : public FuncEval
{
public:
//! Constructor for multiple surfaces.
/*!
* @param k Vector of pointers to InterfaceKinetics objects Each object
* consists of a surface or an edge containing internal degrees of
* freedom representing the concentration of surface adsorbates.
* @param rtol The relative tolerance for the integrator
* @param atol The absolute tolerance for the integrator
* @param maxStepSize The maximum step-size the integrator is allowed to take.
* If zero, this option is disabled.
* @param maxSteps The maximum number of time-steps the integrator can take.
* If not supplied, uses the default value in the CVodesIntegrator (20000).
* @param maxErrTestFails the maximum permissible number of error test failures
* If not supplied, uses the default value in CVODES (7).
*/
ImplicitSurfChem(vector<InterfaceKinetics*> k,
double rtol=1.e-7, double atol=1.e-14,
double maxStepSize=0, size_t maxSteps=20000,
size_t maxErrTestFails=7);
/**
* Must be called before calling method 'advance'
*/
void initialize(double t0=0.0);
/**
* Set the maximum integration step-size. Note, setting this value to zero
* disables this option
*/
void setMaxStepSize(double maxstep=0.0);
/**
* Set the relative and absolute integration tolerances.
*/
void setTolerances(double rtol=1.e-7, double atol=1.e-14);
/**
* Set the maximum number of CVODES integration steps.
*/
void setMaxSteps(size_t maxsteps=20000);
/**
* Set the maximum number of CVODES error test failures
*/
void setMaxErrTestFails(size_t maxErrTestFails=7);
//! Integrate from t0 to t1. The integrator is reinitialized first.
/*!
* This routine does a time accurate solve from t = t0 to t = t1.
* of the surface problem.
*
* @param t0 Initial Time -> this is an input
* @param t1 Final Time -> This is an input
*/
void integrate(double t0, double t1);
//! Integrate from t0 to t1 without reinitializing the integrator.
/*!
* Use when the coverages have not changed from their values on return
* from the last call to integrate or integrate0.
*
* @param t0 Initial Time -> this is an input
* @param t1 Final Time -> This is an input
*/
void integrate0(double t0, double 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 One of the values defined in @ref solvesp_methods.
* The default is -1, which means that the program will decide.
* @param timeScaleOverride When a pseudo 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,
double timeScaleOverride = 1.0);
// overloaded methods of class FuncEval
//! Return the number of equations
size_t neq() const override {
return m_nv;
}
//! Evaluate the value of ydot[k] at the current conditions
/*!
* @param t Time (seconds)
* @param y Vector containing the current solution vector
* @param ydot Output vector containing the value of the
* derivative of the surface coverages.
* @param p Unused parameter pass-through parameter vector
*/
void eval(double t, double* y, double* ydot, double* p) override;
//! Get the current state of the solution vector
/*!
* @param y Value of the solution vector to be used.
* On output, this contains the initial value
* of the solution.
*/
void getState(double* y) override;
/**
* 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(double* 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 double* 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.
*/
vector<InterfaceKinetics*> & getObjects() {
return m_vecKinPtrs;
}
int checkMatch(vector<ThermoPhase*> m_vec, ThermoPhase* thPtr);
void setIOFlag(int ioFlag) {
m_ioFlag = ioFlag;
}
protected:
//! Set the mixture to a state consistent with solution
//! vector y.
/*!
* This function will set the surface site factions in the underlying
* SurfPhase objects to the current value of the solution vector.
*
* @param y Current value of the solution vector. The length is equal to
* the sum of the number of surface sites in all the surface phases.
*/
void updateState(double* y);
//! vector of pointers to surface phases.
vector<SurfPhase*> m_surf;
//! Vector of pointers to bulk phases
vector<ThermoPhase*> m_bulkPhases;
//! vector of pointers to InterfaceKinetics objects
vector<InterfaceKinetics*> m_vecKinPtrs;
//! Vector of number of species in each Surface Phase
vector<size_t> m_nsp;
vector<size_t> m_specStartIndex;
//! Total number of surface species in all surface phases
/*!
* This is the total number of unknowns in m_mode 0 problem
*/
size_t m_nv = 0;
size_t m_numTotalBulkSpecies = 0;
size_t m_numTotalSpecies = 0;
vector<vector<int>> pLocVec;
//! Pointer to the CVODE integrator
unique_ptr<Integrator> m_integ;
double m_atol, m_rtol; // tolerances
double m_maxstep; //!< max step size
size_t m_nmax; //!< maximum number of steps allowed
size_t m_maxErrTestFails; //!< maximum number of error test failures allowed
vector<double> 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<double> m_concSpecies;
vector<double> m_concSpeciesSave;
/**
* Index into the species vector of the kinetics manager,
* pointing to the first species from the surrounding medium.
*/
int m_mediumSpeciesStart = -1;
/**
* Index into the species vector of the kinetics manager, pointing to the
* first species from the condensed phase of the particles.
*/
int m_bulkSpeciesStart = -1;
/**
* Index into the species vector of the kinetics manager, pointing to the
* first species from the surface of the particles
*/
int m_surfSpeciesStart = -1;
/**
* Pointer to the helper method, Placid, which solves the surface problem.
*/
unique_ptr<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 = true;
private:
//! Controls the amount of printing from this routine
//! and underlying routines.
int m_ioFlag = 0;
};
}
#endif
include/cantera/kinetics/solveSP.h
-491
View File
@@ -1,491 +0,0 @@
/**
* @file solveSP.h Header file for implicit surface problem solver (see @ref
* chemkinetics and class @link Cantera::solveSP solveSP@endlink).
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef SOLVESP_H
#define SOLVESP_H
#include "cantera/kinetics/InterfaceKinetics.h"
#include "cantera/numerics/DenseMatrix.h"
//! @defgroup solvesp_methods Surface Problem Solver Methods
//! @ingroup surfSolverGroup
//! @{
//! 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.
const int SFLUX_INITIALIZE = 1;
//! 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.
const int SFLUX_RESIDUAL = 2;
//! 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 because solution variables have
//! been perturbed from nominal values to create Jacobian entries.
const int SFLUX_JACOBIAN = 3;
//! The transient calculation is performed here for an amount of time
//! specified by "time_scale". It is not guaranteed 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.
const int SFLUX_TRANSIENT = 4;
//! @} End of solvesp_method
//! @defgroup solvesp_bulkFunc Surface Problem Bulk Phase Mode
//! Functionality expected from the bulk phase. This changes the equations
//! that will be used to solve for the bulk mole fractions.
//! @ingroup surfSolverGroup
//! @{
//! Deposition of a bulk phase is to be expected. Bulk mole fractions are
//! determined from ratios of growth rates of bulk species.
const int BULK_DEPOSITION = 1;
//! 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_ETCH = 2;
//! @} End of solvesp_bulkFunc
namespace Cantera
{
class ImplicitSurfChem;
//! @addtogroup surfSolverGroup
//! @{
//! 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:
*
* C_i_j | 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 indices 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
*
* Indices 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` parameter,
* that is input to the solution object. This parameter may have one of the
* values defined in @ref solvesp_methods.
*
* ### Pseudo time stepping algorithm:
* The time step is determined from sdot[], 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. See @ref solvesp_bulkFunc. Currently,
* only the default value of BULK_ETCH is supported.
*/
solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc = BULK_ETCH);
//! Destructor
~solveSP() = default;
solveSP(const solveSP&) = delete;
solveSP& operator=(const solveSP&) = delete;
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.
*
* Uses 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.
*
* @param ifunc Determines the type of solution algorithm to be used. See
* @ref solvesp_methods for possible values.
* @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 1 if the surface problem is successfully solved.
* -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 optionally gets called at the start of every
//! invocation
void print_header(int ioflag, int ifunc, double time_scale,
int damping, double reltol, double abstol);
//! 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, size_t iter,
double update_norm, double resid_norm,
bool do_time, bool final=false);
//! 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.
* @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);
/**
* 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 evaluation
/*!
* @param resid output Vector of residuals, length = m_neq
* @param CSolnSP Vector of species concentrations, unknowns in the
* problem, length = m_neq
* @param CSolnOldSP 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 jac Jacobian to be evaluated.
* @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(DenseMatrix& jac, double* resid,
double* CSolnSP,
const double* CSolnSPOld, const bool do_time,
const double deltaT);
//! Vector of interface kinetics objects
/*!
* Each of these is associated with one and only one surface phase.
*/
vector<InterfaceKinetics*> &m_objects;
//! Total number of equations to solve in the implicit problem.
/*!
* Note, this can be zero, and frequently is
*/
size_t m_neq = 0;
//! This variable determines how the bulk phases are to be handled
/*!
* Possible values are given in @ref solvesp_bulkFunc.
*/
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)
*/
size_t m_numSurfPhases = 0;
//! 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.
*/
size_t m_numTotSurfSpecies = 0;
//! 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
*/
vector<size_t> 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
*/
vector<size_t> m_nSpeciesSurfPhase;
//! Vector of surface phase pointers
/*!
* This is created during the constructor
* Length is equal to the number of surface phases, m_numSurfPhases
*/
vector<SurfPhase*> 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.
*/
vector<size_t> m_eqnIndexStartSolnPhase;
//! 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
*/
size_t m_numBulkPhasesSS = 0;
//! Vector of number of species in the m_numBulkPhases phases.
/*!
* Length is number of bulk phases
*/
vector<size_t> m_numBulkSpecies;
//! 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.
*/
size_t m_numTotBulkSpeciesSS = 0;
//! Vector of bulk phase pointers, length is equal to m_numBulkPhases.
vector<ThermoPhase*> 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.
*/
vector<size_t> m_kinSpecIndex;
//! Index between the equation index and the index of the
//! InterfaceKinetics object
/*!
* Length m_neq
*/
vector<size_t> m_kinObjIndex;
//! Vector containing the indices of the largest species
//! in each surface phase
/*!
* `k = m_spSurfLarge[i]` where `k` is the local species index, that is, it
* varies from 0 to (num species in phase - 1) and `i` is the surface
* phase index in the problem. Length is equal to #m_numSurfPhases.
*/
vector<size_t> m_spSurfLarge;
//! Maximum number of species in any single kinetics operator
//! -> also maxed wrt the total # of solution species
size_t m_maxTotSpecies = 0;
//! Temporary vector with length equal to max m_maxTotSpecies
vector<double> m_netProductionRatesSave;
//! Temporary vector with length equal to max m_maxTotSpecies
vector<double> m_numEqn1;
//! Temporary vector with length equal to max m_maxTotSpecies
vector<double> m_numEqn2;
//! Temporary vector with length equal to max m_maxTotSpecies
vector<double> m_CSolnSave;
//! Solution vector. length MAX(1, m_neq)
vector<double> m_CSolnSP;
//! Saved initial solution vector. length MAX(1, m_neq)
vector<double> m_CSolnSPInit;
//! Saved solution vector at the old time step. length MAX(1, m_neq)
vector<double> m_CSolnSPOld;
//! Weights for the residual norm calculation. length MAX(1, m_neq)
vector<double> m_wtResid;
//! Weights for the species concentrations norm calculation
/*!
* length MAX(1, m_neq)
*/
vector<double> 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<double> m_resid;
//! Vector of mole fractions. length m_maxTotSpecies
vector<double> m_XMolKinSpecies;
//! Jacobian. m_neq by m_neq computed Jacobian matrix for the local
//! Newton's method.
DenseMatrix m_Jac;
public:
int m_ioflag = 0;
};
//! @} End of surfSolverGroup
}
#endif
include/cantera/numerics/CVodesIntegrator.h
+1 -2
View File
@@ -20,8 +20,7 @@ namespace Cantera
/**
* Wrapper class for 'cvodes' integrator from LLNL.
*
* See FuncEval.h. Classes that use CVodesIntegrator:
* ImplicitSurfChem, ReactorNet
* See FuncEval.h. Classes that use CVodesIntegrator: ReactorNet
*/
class CVodesIntegrator : public Integrator
{
src/kinetics/ImplicitSurfChem.cpp
-296
View File
@@ -1,296 +0,0 @@
/**
* @file ImplicitSurfChem.cpp
* Definitions for the implicit integration of surface site density equations
* (see @ref kineticsmgr and class
* @link Cantera::ImplicitSurfChem ImplicitSurfChem@endlink).
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#include "cantera/kinetics/ImplicitSurfChem.h"
#include "cantera/kinetics/solveSP.h"
#include "cantera/thermo/SurfPhase.h"
namespace Cantera
{
ImplicitSurfChem::ImplicitSurfChem(
vector<InterfaceKinetics*> k, double rtol, double atol,
double maxStepSize, size_t maxSteps,
size_t maxErrTestFails) :
m_atol(atol),
m_rtol(rtol),
m_maxstep(maxStepSize),
m_nmax(maxSteps),
m_maxErrTestFails(maxErrTestFails)
{
size_t ntmax = 0;
size_t kinSpIndex = 0;
// Loop over the number of surface kinetics objects
for (size_t n = 0; n < k.size(); n++) {
InterfaceKinetics* kinPtr = k[n];
m_vecKinPtrs.push_back(kinPtr);
SurfPhase* surf = dynamic_cast<SurfPhase*>(&k[n]->thermo(0));
if (surf == nullptr) {
throw CanteraError("ImplicitSurfChem::ImplicitSurfChem",
"kinetics manager contains no surface phase");
}
m_surf.push_back(surf);
size_t nsp = m_surf.back()->nSpecies();
m_nsp.push_back(nsp);
m_nv += m_nsp.back();
size_t nt = k[n]->nTotalSpecies();
ntmax = std::max(nt, ntmax);
m_specStartIndex.push_back(kinSpIndex);
kinSpIndex += nsp;
size_t nPhases = kinPtr->nPhases();
vector<int> pLocTmp(nPhases);
pLocTmp[0] = -int(n);
size_t imatch = npos;
for (size_t ip = 1; ip < nPhases; ip++) {
ThermoPhase* thPtr = & kinPtr->thermo(ip);
if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
m_bulkPhases.push_back(thPtr);
nsp = thPtr->nSpecies();
m_numTotalBulkSpecies += nsp;
imatch = m_bulkPhases.size() - 1;
}
pLocTmp[ip] = int(imatch);
}
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.reset(newIntegrator("CVODE"));
// use backward differencing, with a full Jacobian computed
// numerically, and use a Newton linear iterator
m_integ->setMethod(BDF_Method);
m_integ->setLinearSolverType("DENSE");
m_work.resize(ntmax);
}
int ImplicitSurfChem::checkMatch(vector<ThermoPhase*> 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;
}
void ImplicitSurfChem::getState(double* c)
{
size_t loc = 0;
for (size_t n = 0; n < m_surf.size(); n++) {
m_surf[n]->getCoverages(c + loc);
loc += m_nsp[n];
}
}
void ImplicitSurfChem::setMaxStepSize(double maxstep)
{
m_maxstep = maxstep;
if (m_maxstep > 0) {
m_integ->setMaxStepSize(m_maxstep);
}
}
void ImplicitSurfChem::setTolerances(double rtol, double atol)
{
m_rtol = rtol;
m_atol = atol;
m_integ->setTolerances(m_rtol, m_atol);
}
void ImplicitSurfChem::setMaxSteps(size_t maxsteps)
{
m_nmax = maxsteps;
m_integ->setMaxSteps(static_cast<int>(m_nmax));
}
void ImplicitSurfChem::setMaxErrTestFails(size_t maxErrTestFails)
{
m_maxErrTestFails = maxErrTestFails;
m_integ->setMaxErrTestFails(static_cast<int>(m_maxErrTestFails));
}
void ImplicitSurfChem::initialize(double t0)
{
this->setTolerances(m_rtol, m_atol);
this->setMaxStepSize(m_maxstep);
this->setMaxSteps(m_nmax);
this->setMaxErrTestFails(m_maxErrTestFails);
m_integ->initialize(t0, *this);
}
void ImplicitSurfChem::integrate(double t0, double t1)
{
this->initialize(t0);
if (fabs(t1 - t0) < m_maxstep || m_maxstep == 0) {
// limit max step size on this run to t1 - t0
m_integ->setMaxStepSize(t1 - t0);
}
m_integ->integrate(t1);
updateState(m_integ->solution());
}
void ImplicitSurfChem::integrate0(double t0, double t1)
{
m_integ->integrate(t1);
updateState(m_integ->solution());
}
void ImplicitSurfChem::updateState(double* c)
{
size_t loc = 0;
for (size_t n = 0; n < m_surf.size(); n++) {
m_surf[n]->setCoverages(c + loc);
loc += m_nsp[n];
}
}
void ImplicitSurfChem::eval(double time, double* y, double* ydot, double* p)
{
updateState(y); // synchronize the surface state(s) with y
size_t loc = 0;
for (size_t n = 0; n < m_surf.size(); n++) {
double rs0 = 1.0/m_surf[n]->siteDensity();
m_vecKinPtrs[n]->getNetProductionRates(m_work.data());
double sum = 0.0;
for (size_t k = 1; k < m_nsp[n]; k++) {
ydot[k + loc] = m_work[k] * rs0 * m_surf[n]->size(k);
sum -= ydot[k];
}
ydot[loc] = sum;
loc += m_nsp[n];
}
}
void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
double 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 = make_unique<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(m_concSpecies.data());
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 solveSP
// can fail is if concentrations on input are below zero.
bool rset = false;
for (size_t k = 0; k < m_nv; k++) {
if (m_concSpecies[k] < 0.0) {
rset = true;
m_concSpecies[k] = 0.0;
}
}
if (rset) {
setConcSpecies(m_concSpecies.data());
}
m_surfSolver->m_ioflag = m_ioFlag;
// Save the current solution
m_concSpeciesSave = m_concSpecies;
int retn = m_surfSolver->solveSurfProb(ifunc, time_scale, TKelvin, PGas,
reltol, atol);
if (retn != 1) {
// reset the concentrations
m_concSpecies = m_concSpeciesSave;
setConcSpecies(m_concSpeciesSave.data());
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!");
}
}
}
void ImplicitSurfChem::getConcSpecies(double* const vecConcSpecies) const
{
size_t kstart;
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->getConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->getConcentrations(vecConcSpecies + kstart);
kstart += TP_ptr->nSpecies();
}
}
void ImplicitSurfChem::setConcSpecies(const double* const vecConcSpecies)
{
size_t kstart;
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->setConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->setConcentrations(vecConcSpecies + kstart);
kstart += TP_ptr->nSpecies();
}
}
void ImplicitSurfChem::setCommonState_TP(double TKelvin, double PresPa)
{
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}
}
}
src/kinetics/solveSP.cpp
-692
View File
@@ -1,692 +0,0 @@
/**
* @file: solveSP.cpp Implicit surface site concentration solver
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#include "cantera/kinetics/solveSP.h"
#include "cantera/thermo/SurfPhase.h"
#include "cantera/kinetics/ImplicitSurfChem.h"
namespace Cantera
{
// STATIC ROUTINES DEFINED IN THIS FILE
static double calc_damping(double* x, double* dx, size_t dim, int*);
static double calcWeightedNorm(const double [], const double dx[], size_t);
// solveSP Class Definitions
solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
m_objects(surfChemPtr->getObjects()),
m_bulkFunc(bulkFunc)
{
for (size_t n = 0; n < m_objects.size(); n++) {
InterfaceKinetics* kin = m_objects[n];
SurfPhase* sp = dynamic_cast<SurfPhase*>(&kin->thermo(0));
if (sp == nullptr) {
throw CanteraError("solveSP::solveSP",
"InterfaceKinetics object has no surface phase");
}
m_numSurfPhases++;
m_indexKinObjSurfPhase.push_back(n);
m_ptrsSurfPhase.push_back(sp);
size_t nsp = sp->nSpecies();
m_nSpeciesSurfPhase.push_back(nsp);
m_numTotSurfSpecies += nsp;
}
if (bulkFunc == BULK_DEPOSITION) {
m_neq = m_numTotSurfSpecies + m_numTotBulkSpeciesSS;
} else {
m_neq = m_numTotSurfSpecies;
}
for (size_t n = 0; n < m_numSurfPhases; n++) {
size_t tsp = m_objects[n]->nTotalSpecies();
m_maxTotSpecies = std::max(m_maxTotSpecies, tsp);
}
m_maxTotSpecies = std::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);
size_t kindexSP = 0;
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
size_t iKinObject = m_indexKinObjSurfPhase[isp];
InterfaceKinetics* kin = m_objects[iKinObject];
size_t kstart = kin->kineticsSpeciesIndex(0, 0);
size_t nsp = m_nSpeciesSurfPhase[isp];
m_eqnIndexStartSolnPhase[isp] = kindexSP;
for (size_t k = 0; k < nsp; k++, kindexSP++) {
m_kinSpecIndex[kindexSP] = kstart + k;
m_kinObjIndex[kindexSP] = isp;
}
}
// Dimension solution vector
size_t dim1 = std::max<size_t>(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_Jac.resize(dim1, dim1, 0.0);
}
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 label_t=-1; // Species IDs for time control
int label_d = -1; // Species IDs for damping control
int label_t_old=-1;
double label_factor = 1.0;
int iter=0; // iteration number on nonlinear solver
int iter_max=1000; // maximum number of nonlinear iterations
double deltaT = 1.0E-10; // Delta time step
double damp=1.0;
double inv_t = 0.0;
double t_real = 0.0, update_norm = 1.0E6;
bool do_time = false, not_converged = true;
m_ioflag = std::min(m_ioflag, 1);
// 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.
size_t loc = 0;
for (size_t n = 0; n < m_numSurfPhases; n++) {
m_ptrsSurfPhase[n]->getConcentrations(m_numEqn1.data());
for (size_t k = 0; k < m_nSpeciesSurfPhase[n]; k++) {
m_CSolnSP[loc] = m_numEqn1[k];
loc++;
}
}
m_CSolnSPInit = m_CSolnSP;
// Calculate the largest species in each phase
evalSurfLarge(m_CSolnSP.data());
if (m_ioflag) {
print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol);
}
// 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
m_CSolnSPOld = m_CSolnSP;
// Evaluate the largest surface species for each surface phase every
// 5 iterations.
if (iter%5 == 4) {
evalSurfLarge(m_CSolnSP.data());
}
// 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;
}
double tmp = calc_t(m_netProductionRatesSave.data(),
m_XMolKinSpecies.data(),
&label_t, &label_t_old, &label_factor, m_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_Jac, m_resid.data(), m_CSolnSP.data(),
m_CSolnSPOld.data(), do_time, deltaT);
// Calculate the weights. Make sure the calculation is carried out on
// the first iteration.
if (iter%4 == 1) {
calcWeights(m_wtSpecies.data(), m_wtResid.data(),
m_Jac, m_CSolnSP.data(), abstol, reltol);
}
// Find the weighted norm of the residual
double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
// Solve Linear system. The solution is in m_resid
solve(m_Jac, m_resid.data());
// 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.
damp = calc_damping(m_CSolnSP.data(), m_resid.data(), m_neq, &label_d);
// Calculate the weighted norm of the update vector Here, resid is the
// delta of the solution, in concentration units.
update_norm = calcWeightedNorm(m_wtSpecies.data(),
m_resid.data(), m_neq);
// Update the solution vector and real time Crop the concentrations to
// zero.
for (size_t irow = 0; irow < m_neq; irow++) {
m_CSolnSP[irow] -= damp * m_resid[irow];
}
for (size_t irow = 0; irow < m_neq; irow++) {
m_CSolnSP[irow] = std::max(0.0, m_CSolnSP[irow]);
}
updateState(m_CSolnSP.data());
if (do_time) {
t_real += damp/inv_t;
}
if (m_ioflag) {
printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
update_norm, resid_norm, 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;
}
} 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 && m_ioflag) {
writelog("#$#$#$# Error in solveSP $#$#$#$ \n");
writelogf("Newton iter on surface species did not converge, "
"update_norm = %e \n", update_norm);
writelog("Continuing anyway\n");
}
// Decide on what to return in the solution vector. Right now, will always
// return the last solution no matter how bad
if (m_ioflag) {
fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(),
false, deltaT);
double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
update_norm, resid_norm, do_time, true);
}
// Return with the appropriate flag
if (update_norm > 1.0) {
return -1;
}
return 1;
}
void solveSP::updateState(const double* CSolnSP)
{
vector<double> X;
size_t loc = 0;
for (size_t n = 0; n < m_numSurfPhases; n++) {
X.resize(m_nSpeciesSurfPhase[n]);
for (size_t k = 0; k < X.size(); k++) {
X[k] = CSolnSP[loc + k] / m_ptrsSurfPhase[n]->siteDensity();
}
m_ptrsSurfPhase[n]->setMoleFractions_NoNorm(X.data());
loc += m_nSpeciesSurfPhase[n];
}
}
void solveSP::updateMFSolnSP(double* XMolSolnSP)
{
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
size_t keqnStart = m_eqnIndexStartSolnPhase[isp];
m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart);
}
}
void solveSP::updateMFKinSpecies(double* XMolKinSpecies, int isp)
{
InterfaceKinetics* kin = m_objects[isp];
for (size_t iph = 0; iph < kin->nPhases(); iph++) {
size_t ksi = kin->kineticsSpeciesIndex(0, iph);
kin->thermo(iph).getMoleFractions(XMolKinSpecies + ksi);
}
}
void solveSP::evalSurfLarge(const double* CSolnSP)
{
size_t kindexSP = 0;
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
double Clarge = CSolnSP[kindexSP];
m_spSurfLarge[isp] = 0;
kindexSP++;
for (size_t k = 1; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
if (CSolnSP[kindexSP] > Clarge) {
Clarge = CSolnSP[kindexSP];
m_spSurfLarge[isp] = k;
}
}
}
}
void solveSP::fun_eval(double* resid, const double* CSoln, const double* CSolnOld,
const bool do_time, const double deltaT)
{
size_t k;
double lenScale = 1.0E-9;
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) {
size_t kindexSP = 0;
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
size_t nsp = m_nSpeciesSurfPhase[isp];
InterfaceKinetics* kinPtr = m_objects[isp];
size_t kins = kindexSP;
kinPtr->getNetProductionRates(m_netProductionRatesSave.data());
for (k = 0; k < nsp; k++, kindexSP++) {
resid[kindexSP] =
(CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
- m_netProductionRatesSave[k];
}
size_t kspecial = kins + m_spSurfLarge[isp];
double sd = m_ptrsSurfPhase[isp]->siteDensity();
resid[kspecial] = sd;
for (k = 0; k < nsp; k++) {
resid[kspecial] -= CSoln[kins + k];
}
}
} else {
size_t kindexSP = 0;
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
size_t nsp = m_nSpeciesSurfPhase[isp];
InterfaceKinetics* kinPtr = m_objects[isp];
size_t kins = kindexSP;
kinPtr->getNetProductionRates(m_netProductionRatesSave.data());
for (k = 0; k < nsp; k++, kindexSP++) {
resid[kindexSP] = - m_netProductionRatesSave[k];
}
size_t kspecial = kins + m_spSurfLarge[isp];
double sd = m_ptrsSurfPhase[isp]->siteDensity();
resid[kspecial] = sd;
for (k = 0; k < nsp; k++) {
resid[kspecial] -= CSoln[kins + k];
}
}
}
if (m_bulkFunc == BULK_DEPOSITION) {
size_t kindexSP = m_numTotSurfSpecies;
for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
double* XBlk = m_numEqn1.data();
size_t nsp = m_nSpeciesSurfPhase[isp];
double grRate = 0.0;
for (k = 0; k < nsp; k++) {
if (m_netProductionRatesSave[k] > 0.0) {
grRate += m_netProductionRatesSave[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[k] > 0.0) {
resid[kindexSP + k] = XBlk[k] * grRate
- m_netProductionRatesSave[k];
} else {
resid[kindexSP + k] = XBlk[k] * grRate;
}
}
} else {
grRate = 1.0E-6;
//! @todo the appearance of k in this formula is suspicious
grRate += fabs(m_netProductionRatesSave[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;
}
}
}
}
void solveSP::resjac_eval(DenseMatrix& jac, double resid[], double CSoln[],
const double CSolnOld[], const bool do_time,
const double deltaT)
{
size_t kColIndex = 0;
// Calculate the residual
fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
// Now we will look over the columns perturbing each unknown.
for (size_t jsp = 0; jsp < m_numSurfPhases; jsp++) {
size_t nsp = m_nSpeciesSurfPhase[jsp];
double sd = m_ptrsSurfPhase[jsp]->siteDensity();
for (size_t kCol = 0; kCol < nsp; kCol++) {
double cSave = CSoln[kColIndex];
double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
CSoln[kColIndex] += dc;
fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
for (size_t i = 0; i < m_neq; i++) {
jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
}
CSoln[kColIndex] = cSave;
kColIndex++;
}
}
if (m_bulkFunc == BULK_DEPOSITION) {
for (size_t jsp = 0; jsp < m_numBulkPhasesSS; jsp++) {
size_t nsp = m_numBulkSpecies[jsp];
double sd = m_bulkPhasePtrs[jsp]->molarDensity();
for (size_t kCol = 0; kCol < nsp; kCol++) {
double cSave = CSoln[kColIndex];
double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
CSoln[kColIndex] += dc;
fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
for (size_t i = 0; i < m_neq; i++) {
jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
}
CSoln[kColIndex] = cSave;
kColIndex++;
}
}
}
}
/**
* 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.
*/
static double calc_damping(double x[], double dxneg[], size_t dim, int* label)
{
const double APPROACH = 0.80;
double damp = 1.0;
static double damp_old = 1.0; //! @todo this variable breaks thread safety
*label = -1;
for (size_t i = 0; i < dim; i++) {
// Calculate the new suggested new value of x[i]
double 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
double xtop = 1.0 - 0.1*fabs(1.0-x[i]);
double xbot = fabs(x[i]*0.1) - 1.0e-16;
if (xnew > xtop) {
damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
*label = int(i);
} else if (xnew < xbot) {
damp = APPROACH * x[i] / dxneg[i];
*label = int(i);
} else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
*label = int(i);
}
}
damp = std::max(damp, 1e-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 */
/**
* 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[], size_t dim)
{
double norm = 0.0;
if (dim == 0) {
return 0.0;
}
for (size_t i = 0; i < dim; i++) {
norm += pow(dx[i] / wtX[i], 2);
}
return sqrt(norm/dim);
}
void solveSP::calcWeights(double wtSpecies[], double wtResid[],
const Array2D& Jac, const double CSoln[],
const double abstol, const double reltol)
{
// First calculate the weighting factor for the concentrations of the
// surface species and bulk species.
size_t kindex = 0;
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
double sd = m_ptrsSurfPhase[isp]->siteDensity();
for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindex++) {
wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
}
}
if (m_bulkFunc == BULK_DEPOSITION) {
for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
double sd = m_bulkPhasePtrs[isp]->molarDensity();
for (size_t k = 0; k < m_numBulkSpecies[isp]; 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 (size_t k = 0; k < m_neq; k++) {
wtResid[k] = 0.0;
for (size_t jcol = 0; jcol < m_neq; jcol++) {
wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
}
}
}
double solveSP::calc_t(double netProdRateSolnSP[], double XMolSolnSP[], int* label,
int* label_old, double* label_factor, int ioflag)
{
double inv_timeScale = 1.0E-10;
size_t kindexSP = 0;
*label = 0;
updateMFSolnSP(XMolSolnSP);
for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
// Get the interface kinetics associated with this surface
InterfaceKinetics* kin = m_objects[isp];
kin->getNetProductionRates(m_numEqn1.data());
double sden = kin->thermo(0).molarDensity();
for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
netProdRateSolnSP[kindexSP] = m_numEqn1[k];
double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10);
tmp *= sden;
tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
if (netProdRateSolnSP[kindexSP]> 0.0) {
tmp /= 100.;
}
if (tmp > inv_timeScale) {
inv_timeScale = tmp;
*label = int(kindexSP);
}
}
}
// 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;
}
return inv_timeScale / *label_factor;
} // calc_t
void solveSP::print_header(int ioflag, int ifunc, double time_scale,
int damping, double reltol, double abstol)
{
if (ioflag) {
writelog("\n================================ SOLVESP CALL SETUP "
"========================================\n");
if (ifunc == SFLUX_INITIALIZE) {
writelog("\n SOLVESP Called with Initialization turned on\n");
writelogf(" Time scale input = %9.3e\n", time_scale);
} else if (ifunc == SFLUX_RESIDUAL) {
writelog("\n SOLVESP Called to calculate steady state residual\n");
writelog(" from a good initial guess\n");
} else if (ifunc == SFLUX_JACOBIAN) {
writelog("\n SOLVESP Called to calculate steady state Jacobian\n");
writelog(" from a good initial guess\n");
} else if (ifunc == SFLUX_TRANSIENT) {
writelog("\n SOLVESP Called to integrate surface in time\n");
writelogf(" for a total of %9.3e sec\n", time_scale);
} else {
throw CanteraError("solveSP::print_header",
"Unknown ifunc flag = {}", ifunc);
}
if (m_bulkFunc == BULK_DEPOSITION) {
writelog(" The composition of the Bulk Phases will be calculated\n");
} else if (m_bulkFunc == BULK_ETCH) {
writelog(" Bulk Phases have fixed compositions\n");
} else {
throw CanteraError("solveSP::print_header",
"Unknown bulkFunc flag = {}", m_bulkFunc);
}
if (damping) {
writelog(" Damping is ON \n");
} else {
writelog(" Damping is OFF \n");
}
writelogf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
}
if (ioflag == 1) {
writelog("\n\n\t Iter Time Del_t Damp DelX "
" Resid Name-Time Name-Damp\n");
writelog("\t -----------------------------------------------"
"------------------------------------\n");
}
}
void solveSP::printIteration(int ioflag, double damp, int label_d,
int label_t, double inv_t, double t_real,
size_t iter, double update_norm,
double resid_norm, bool do_time, bool final)
{
if (ioflag == 1) {
if (final) {
writelogf("\tFIN%3d ", iter);
} else {
writelogf("\t%6d ", iter);
}
if (do_time) {
writelogf("%9.4e %9.4e ", t_real, 1.0/inv_t);
} else {
writeline(' ', 22, false);
}
if (damp < 1.0) {
writelogf("%9.4e ", damp);
} else {
writeline(' ', 11, false);
}
writelogf("%9.4e %9.4e", update_norm, resid_norm);
if (do_time) {
size_t k = m_kinSpecIndex[label_t];
size_t isp = m_kinObjIndex[label_t];
writelog(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
} else {
writeline(' ', 16, false);
}
if (label_d >= 0) {
size_t k = m_kinSpecIndex[label_d];
size_t isp = m_kinObjIndex[label_d];
writelogf(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
}
if (final) {
writelog(" -- success");
}
writelog("\n");
}
} // printIteration
}
test_problems/.gitignore
-1
View File
@@ -31,5 +31,4 @@ diamondSurf/runDiamond_output.out
dustyGasTransport/dustyGasTransport
pureFluidTest/pureFluid
surfSolverTest/surfaceSolver
surfSolverTest/surfaceSolver2
stoichSolidKinetics/stoichSolidKinetics
test_problems/SConscript
-4
View File
@@ -194,10 +194,6 @@ CompileAndTest('surfSolver', 'surfSolverTest', 'surfaceSolver', None,
comparisons=[('results_blessed.txt', 'results.txt')],
artifacts=['results.txt'],
extensions=['^surfaceSolver.cpp'])
CompileAndTest('surfSolver2', 'surfSolverTest', 'surfaceSolver2', None,
comparisons=[('results2_blessed.txt', 'results2.txt')],
artifacts=['results2.txt'],
extensions=['^surfaceSolver2.cpp'])
CompileAndTest('VCS-NaCl', 'VCSnonideal/NaCl_equil',
'nacl_equil', 'good_out.txt',
options='-d 3',
test_problems/surfSolverTest/results2_blessed.txt
-44
View File
@@ -1,44 +0,0 @@
Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05
Gas Phase: gas (2)
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 H 8.69602e-05 0.00999001 -2.365619e-03
1 N2 0.00773858 0.889011 0.000000e+00
2 CH3 0 0 0.000000e+00
3 CH4 8.69601e-05 0.00999001 0.000000e+00
4 C2H2 8.69602e-05 0.00999001 -1.007791e-04
5 H2 0.000347841 0.03996 1.282498e-03
6 OH 8.69602e-07 9.99001e-05 -7.542502e-05
7 H2O 0.000347841 0.03996 3.880308e-05
8 CO 0 0 3.843405e-05
9 O2 8.69602e-06 0.000999001 -9.060640e-07
Sum of gas mole fractions= 1
Bulk Phase: soot (11)
Bulk Temperature = 1.4e+03
Bulk Pressure = 1.01e+05
Name Conc MoleF SrcRate
(kmol/m^3) (kmol/m^2/s)
0 CB-CB3 293.065 1 1.631241e-04
Bulk Weight Growth Rate = 0.00196 kg/m^2/s
Bulk Growth Rate = 5.57e-07 m/s
Bulk Growth Rate = 2e+03 microns / hour
Density of bulk phase = 3.52e+03 kg / m^3
= 3.52 gm / cm^3
Sum of bulk mole fractions= 1
Surface Phase: soot_interface (0)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 0.000000e+00
1 Csoot-H 0.981532 0.000000e+00
Sum of coverages = 1
Surface Phase: soot_interface (0)
Surface Temperature = 1.4e+03
Surface Pressure = 1.01e+05
Name Coverage SrcRate
0 Csoot-* 0.0184677 0.000000e+00
1 Csoot-H 0.981532 0.000000e+00
Sum of coverages = 1
test_problems/surfSolverTest/surfaceSolver2.cpp
-279
View File
@@ -1,279 +0,0 @@
/**
* @file surfaceSolver2.cpp
*
*/
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#define MSSIZE 200
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/kinetics/KineticsFactory.h"
#include "cantera/kinetics/ImplicitSurfChem.h"
#include <iostream>
#include <fstream>
using namespace std;
using namespace Cantera;
void printGas(ostream& oooo, shared_ptr<ThermoPhase> gasTP,
shared_ptr<InterfaceKinetics> iKin_ptr, double* src)
{
double x[MSSIZE];
double C[MSSIZE];
oooo.precision(3);
string gasPhaseName = "gas";
gasTP->getMoleFractions(x);
gasTP->getConcentrations(C);
double Temp = gasTP->temperature();
double p = gasTP->pressure();
oooo << "Gas Temperature = " << Temp << endl;
oooo << "Gas Pressure = " << p << endl;
size_t nGas = iKin_ptr->phaseIndex(gasPhaseName);
size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nGas);
oooo << "Gas Phase: " << gasPhaseName << " "
<< "(" << kstart << ")" << endl;
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
size_t nspGas = gasTP->nSpecies();
for (size_t k = 0; k < nspGas; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, nGas);
fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14g} {:14e}\n",
k, gasTP->speciesName(k), C[k], x[k], src[kstart]);
sum += x[k];
}
oooo << "Sum of gas mole fractions= " << sum << endl;
oooo << endl;
}
void printBulk(ostream& oooo, shared_ptr<ThermoPhase> bulkPhaseTP,
shared_ptr<InterfaceKinetics> iKin_ptr, double* src)
{
double x[MSSIZE];
double C[MSSIZE];
oooo.precision(3);
string bulkParticlePhaseName = bulkPhaseTP->name();
bulkPhaseTP->getMoleFractions(x);
bulkPhaseTP->getConcentrations(C);
size_t nBulk = iKin_ptr->phaseIndex(bulkParticlePhaseName);
size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nBulk);
double dens = bulkPhaseTP->density();
oooo << "Bulk Phase: " << bulkParticlePhaseName << " "
<< "(" << kstart << ")" << endl;
double Temp = bulkPhaseTP->temperature();
double p = bulkPhaseTP->pressure();
oooo << "Bulk Temperature = " << Temp << endl;
oooo << "Bulk Pressure = " << p << endl;
oooo << " Name "
<< " Conc MoleF SrcRate " << endl;
oooo << " "
<< " (kmol/m^3) (kmol/m^2/s) " << endl;
double sum = 0.0;
double Wsum = 0.0;
const vector<double>& molecW = bulkPhaseTP->molecularWeights();
size_t nspBulk = bulkPhaseTP->nSpecies();
for (size_t k = 0; k < nspBulk; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, nBulk);
fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14g} {:14e}\n",
k, bulkPhaseTP->speciesName(k), C[k], x[k], src[kstart]);
sum += x[k];
Wsum += src[kstart] * molecW[k];
}
oooo << "Bulk Weight Growth Rate = " << Wsum << " kg/m^2/s" << endl;
double gr = Wsum / dens;
oooo << "Bulk Growth Rate = " << gr << " m/s" << endl;
oooo << "Bulk Growth Rate = " << gr * 1.0E6 * 3600.
<< " microns / hour" << endl;
oooo << "Density of bulk phase = " << dens << " kg / m^3 "<< endl;
oooo << " = " << dens / 1.0E3
<<" gm / cm^3 " << endl;
oooo << "Sum of bulk mole fractions= " << sum << endl;
oooo << endl;
}
void printSurf(ostream& oooo, shared_ptr<ThermoPhase> surfPhaseTP,
shared_ptr<InterfaceKinetics> iKin_ptr, double* src)
{
double x[MSSIZE];
oooo.precision(3);
string surfParticlePhaseName = surfPhaseTP->name();
surfPhaseTP->getMoleFractions(x);
size_t nSurf = iKin_ptr->phaseIndex(surfParticlePhaseName);
size_t kstart = iKin_ptr->kineticsSpeciesIndex(0, nSurf);
oooo << "Surface Phase: " << surfParticlePhaseName
<< " (" << kstart << ")" << endl;
double Temp = surfPhaseTP->temperature();
double p = surfPhaseTP->pressure();
oooo << "Surface Temperature = " << Temp << endl;
oooo << "Surface Pressure = " << p << endl;
oooo << " Name "
<< " Coverage SrcRate " << endl;
double sum = 0.0;
size_t nspSurf = surfPhaseTP->nSpecies();
for (size_t k = 0; k < nspSurf; k++) {
kstart = iKin_ptr->kineticsSpeciesIndex(k, nSurf);
double srcK = src[kstart];
if (fabs(srcK) < 1.0E-8) {
srcK = 0.0;
}
fmt::print(oooo, "{:4d} {:>24s} {:14g} {:14e}\n",
k, surfPhaseTP->speciesName(k), x[k], srcK);
sum += x[k];
}
oooo << "Sum of coverages = " << sum << endl;
}
int main(int argc, char** argv)
{
string infile = "haca2.yaml";
string gasPhaseName = "gas";
string bulkParticlePhaseName = "soot";
string surfParticlePhaseName = "soot_interface";
int ioflag = 1;
try {
auto gasTP = newThermo(infile, gasPhaseName);
size_t nspGas = gasTP->nSpecies();
cout << "Number of species = " << nspGas << endl;
auto bulkPhaseTP = newThermo(infile, bulkParticlePhaseName);
size_t nspBulk = bulkPhaseTP->nSpecies();
cout << "Number of species in bulk phase named " <<
bulkParticlePhaseName << " = " << nspBulk << endl;
auto surfPhaseTP = newThermo(infile, surfParticlePhaseName);
size_t nsp_d100 = surfPhaseTP->nSpecies();
cout << "Number of species in surface phase, " << surfParticlePhaseName
<< " = " << nsp_d100 << endl;
auto kin = newKinetics({surfPhaseTP, gasTP, bulkPhaseTP}, infile);
auto iKin_ptr = dynamic_pointer_cast<InterfaceKinetics>(kin);
size_t nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
ofstream ofile("results2.txt");
// create a second copy of the same surface phase
// (this is a made up problem btw to check the software capability)
auto surfPhaseTP2 = newThermo(infile, surfParticlePhaseName);
size_t nsp2 = surfPhaseTP2->nSpecies();
string pname = surfPhaseTP2->name();
cout << "Number of species in 2nd surface phase, " << pname
<< " = " << nsp2 << endl;
// create the second InterfaceKinetics object based on the
// second surface phase.
auto kin2 = newKinetics({surfPhaseTP2, gasTP, bulkPhaseTP}, infile);
auto iKin2_ptr = dynamic_pointer_cast<InterfaceKinetics>(kin2);
nr = iKin_ptr->nReactions();
cout << "Number of reactions = " << nr << endl;
double x[MSSIZE];
/*
* Set-up the Surface Problem
* This problem will consist of 2 identical InterfaceKinetics objects
*/
vector<InterfaceKinetics*> vecKinPtrs { iKin_ptr.get(), iKin2_ptr.get() };
// Create the ImplicitSurfChem problem
// Initialize it and call the pseudo steadystate capability.
ImplicitSurfChem* surfaceProb = new ImplicitSurfChem(vecKinPtrs);
surfaceProb->initialize();
surfaceProb->solvePseudoSteadyStateProblem(ioflag);
/*
* Download the source terms for the rate equations
*/
double src[MSSIZE];
double src2[MSSIZE];
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
printGas(ofile, gasTP, iKin_ptr, src);
printBulk(ofile, bulkPhaseTP, iKin_ptr, src);
printSurf(ofile, surfPhaseTP, iKin_ptr, src) ;
printSurf(ofile, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
/*
* Set the Gas State:
* -> note that the states are set in the input file too
*/
double pres = gasTP->pressure();
gasTP->getMoleFractions(x);
double tmp = 0.3 * std::min(x[0], x[1]);
x[0] += tmp;
x[1] -= tmp;
gasTP->setMoleFractions(x);
gasTP->setPressure(pres);
surfaceProb->solvePseudoSteadyStateProblem(ioflag);
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Tweak the inputs and do a quick calculation */
/****************************************************************************/
pres = surfPhaseTP->pressure();
double temp = surfPhaseTP->temperature();
temp += 95;
surfPhaseTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem(ioflag);
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
/*****************************************************************************/
/* Now Don't Tweak the inputs at all */
/****************************************************************************/
gasTP->setState_TP(temp, pres);
surfaceProb->solvePseudoSteadyStateProblem(ioflag);
iKin_ptr->getNetProductionRates(src);
iKin2_ptr->getNetProductionRates(src2);
printGas(cout, gasTP, iKin_ptr, src);
printBulk(cout, bulkPhaseTP, iKin_ptr, src);
printSurf(cout, surfPhaseTP, iKin_ptr, src) ;
printSurf(cout, surfPhaseTP2, iKin2_ptr, src2) ;
delete surfaceProb;
surfaceProb = 0;
iKin_ptr = 0;
iKin2_ptr = 0;
appdelete();
} catch (CanteraError& err) {
std::cout << err.what() << std::endl;
}
return 0;
}
/***********************************************************/