From 2cc7c139b67957af6b8f9b2d678fc90e38182b5a Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 24 Dec 2007 15:27:37 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/src/kinetics/EdgeKinetics.h | 10 +- Cantera/src/kinetics/GRI_30_Kinetics.cpp | 1 + Cantera/src/kinetics/GRI_30_Kinetics.h | 1 + Cantera/src/kinetics/InterfaceKinetics.h | 9 + Cantera/src/kinetics/Kinetics.h | 5 + Cantera/src/kinetics/solveSP.cpp | 2495 +++++++++++----------- Cantera/src/spectra/LineBroadener.cpp | 117 +- Cantera/src/spectra/LineBroadener.h | 54 +- Cantera/src/spectra/rotor.cpp | 21 +- Cantera/src/spectra/rotor.h | 31 +- Cantera/src/thermo/AdsorbateThermo.h | 173 +- Cantera/src/thermo/Phase.h | 48 +- Cantera/src/thermo/PureFluidPhase.h | 20 +- Cantera/src/thermo/SpeciesThermo.h | 2 +- Cantera/src/thermo/ThermoPhase.h | 5 +- Cantera/src/transport/TransportBase.h | 8 + 16 files changed, 1588 insertions(+), 1412 deletions(-) diff --git a/Cantera/src/kinetics/EdgeKinetics.h b/Cantera/src/kinetics/EdgeKinetics.h index f1eaf2de8..00bb6d15f 100644 --- a/Cantera/src/kinetics/EdgeKinetics.h +++ b/Cantera/src/kinetics/EdgeKinetics.h @@ -1,7 +1,11 @@ /** * @file EdgeKinetics.h * - * $Author$ + * @ingroup chemkinetics + * @ingroup electrochem + */ + +/* $Author$ * $Revision$ * $Date$ */ @@ -16,6 +20,10 @@ namespace Cantera { + /** + * Heterogeneous reactions at one-dimensional interfaces between + * multiple adjacent two-dimensional surfaces. + */ class EdgeKinetics : public InterfaceKinetics { public: diff --git a/Cantera/src/kinetics/GRI_30_Kinetics.cpp b/Cantera/src/kinetics/GRI_30_Kinetics.cpp index f71429275..7bd3cbe74 100755 --- a/Cantera/src/kinetics/GRI_30_Kinetics.cpp +++ b/Cantera/src/kinetics/GRI_30_Kinetics.cpp @@ -1,6 +1,7 @@ /** * @file GRI_30_Kinetics.cpp * + * @ingroup chemkinetics */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/GRI_30_Kinetics.h b/Cantera/src/kinetics/GRI_30_Kinetics.h index c9dd601d8..0d6d16fca 100755 --- a/Cantera/src/kinetics/GRI_30_Kinetics.h +++ b/Cantera/src/kinetics/GRI_30_Kinetics.h @@ -2,6 +2,7 @@ * * @file GRI_30_Kinetics.h * + * @ingroup chemkinetics */ // Copyright 2001 California Institute of Technology diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index e6bc65146..3dd7610f2 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -1,6 +1,7 @@ /** * @file InterfaceKinetics.h * + * @ingroup chemkinetics */ /* * $Author$ @@ -10,6 +11,13 @@ // Copyright 2001 California Institute of Technology +/** + * @defgroup electrochem Electrochemistry + * + * Support for electrochemical reaction kinetics. + * + * @ingroup chemkinetics + */ #ifndef CT_IFACEKINETICS_H #define CT_IFACEKINETICS_H @@ -67,6 +75,7 @@ namespace Cantera { /// reactions are assumed to occur at a 2D interface between two /// 3D phases. /// + /// @ingroup chemkinetics class InterfaceKinetics : public Kinetics { public: diff --git a/Cantera/src/kinetics/Kinetics.h b/Cantera/src/kinetics/Kinetics.h index fe2651104..de1ffb924 100755 --- a/Cantera/src/kinetics/Kinetics.h +++ b/Cantera/src/kinetics/Kinetics.h @@ -23,6 +23,10 @@ namespace Cantera { // forward references class ReactionData; + /** + * @defgroup chemkinetics Chemical Kinetics + */ + /// @defgroup kineticsmgr Kinetics Managers /// @section kinmodman Models and Managers /// @@ -106,6 +110,7 @@ namespace Cantera { /// in the reactions), the next 3 will be for phase 'b', and finally the /// net production rates for the surface species will occupy the last /// 5 locations. + /// @ingroup chemkinetics //! Public interface for kinetics managers. diff --git a/Cantera/src/kinetics/solveSP.cpp b/Cantera/src/kinetics/solveSP.cpp index a70dea624..0d0985b0b 100644 --- a/Cantera/src/kinetics/solveSP.cpp +++ b/Cantera/src/kinetics/solveSP.cpp @@ -26,25 +26,25 @@ using namespace std; namespace Cantera { - /*************************************************************************** - * STATIC ROUTINES DEFINED IN THIS FILE - ***************************************************************************/ + /*************************************************************************** + * 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); + 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 - *****************************************************************************/ + /*************************************************************************** + * 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 */ @@ -54,1392 +54,1389 @@ namespace Cantera { # define DAMPING true #endif - /*************************************************************************** - * solveSP Class Definitinos - ***************************************************************************/ + /*************************************************************************** + * 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), - m_ioflag(0) - { + // 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), + m_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_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; + 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; - //} + } + /* + * 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]; - //} - //} - } + //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; - } + 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_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_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_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); + 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; - //} - //} - } + 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); + // 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); + 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; + // Empty destructor + solveSP::~solveSP() { } - 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; + /* + * 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; + double t1; #endif #else - if (m_ioflag > 1) { - m_ioflag = 1; - } + if (m_ioflag > 1) { + m_ioflag = 1; + } #endif #ifdef DEBUG_SOLVESP #ifdef DEBUG_SOLVESP_TIME - Cantera::clockWC wc; - if (m_ioflag) t1 = wc.secondsWC(); + Cantera::clockWC wc; + if (m_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; + /* + * 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_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()); + if (m_bulkFunc == BULK_DEPOSITION) { + //for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + //ThermoPhase *bf_ptr = m_bulkPhasePtrs[isp]; + //bf_ptr->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)); + // 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 (m_ioflag) { - print_header(m_ioflag, ifunc, time_scale, DAMPING, reltol, abstol, - TKelvin, PGas, DATA_PTR(m_netProductionRatesSave), - DATA_PTR(m_XMolKinSpecies)); - } + if (m_ioflag) { + print_header(m_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; - } + /* + * 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()); + /* ------------------------------------------------------------------ + * 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)); - } + /* + * 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, 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; - } + /* + * 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, 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 - */ + /* + * 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; + 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); + /* + * 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), + /* + * 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); + /* + * Find the weighted norm of the residual + */ + resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid), + DATA_PTR(m_resid), m_neq); #ifdef DEBUG_SOLVESP - if (m_ioflag > 1) { - printIterationHeader(m_ioflag, damp, inv_t, t_real, iter, do_time); - /* - * Print out the residual and jacobian - */ - printResJac(m_ioflag, m_neq, m_Jac, DATA_PTR(m_resid), + if (m_ioflag > 1) { + printIterationHeader(m_ioflag, damp, inv_t, t_real, iter, do_time); + /* + * Print out the residual and jacobian + */ + printResJac(m_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 (m_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; + /* + * 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) { + 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 (m_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 (m_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; + /* print out some helpful info */ + if (m_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 (m_ioflag) { - printf("\nResidual is small, forcing convergence!\n"); - } + if (m_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. - */ + /* + * 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); + 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)); + /* + * 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 (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, - 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 (m_ioflag) { + printIteration(m_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; + 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 (m_ioflag > 1) { - printf("\t\tSwitching to steady solve.\n"); - } + if (m_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 */ + } + } + 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 (m_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"); - } - } + /* + * End Newton's method. If not converged, print error message and + * recalculate sdot's at equal site fractions. + */ + if (not_converged) { + if (m_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 (m_ioflag) { - printf("\nEnd of solve, time used: %e\n", wc.secondsWC()-t1); - } + if (m_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 (m_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(m_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); - } + /* + * 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(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(m_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 with the appropriate flag + */ + if (update_norm > 1.0) { + return -1; + } + 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]; + /* + * 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]; + // } + //} } - //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); + /* + * 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); + // } + //} } - //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 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; - } - } + /* + * 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. - */ + /* + * 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]; - } + 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]; - } - } - } + 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; - } - } + 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 + * Calculate the Jacobian and residual + * + * @internal + * This routine uses the m_numEqn2 vector + * as temporary internal storage. */ - 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++; - } - } + 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++; - } - } + 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) + 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. - */ + /* 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; + { + int i; + double damp = 1.0, xnew, xtop, xbot; + static double damp_old = 1.0; - *label = -1; + *label = -1; - for (i = 0; i < dim; i++) { + for (i = 0; i < dim; i++) { - /* - * Calculate the new suggested new value of x[i] - */ + /* + * Calculate the new suggested new value of x[i] + */ - xnew = x[i] - damp * dxneg[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 - */ + /* + * 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; - } - } + 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 < 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; - } + if (damp > damp_old*3) { + damp = damp_old*3; + *label = -1; + } - /* - * Save old value of the damping parameter for use - * in subsequent calls. - */ + /* + * Save old value of the damping parameter for use + * in subsequent calls. + */ - damp_old = damp; - return damp; + damp_old = damp; + return damp; - } /* calc_damping */ + } /* 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. + * This function calculates the norm of an update, dx[], + * based on the weighted values of x. */ - 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]); - } - } - } + 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)); + } - /* - * 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]; + /* + * 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]; + // 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); + // 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)); + 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; - } - } - } + 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; + /* + * 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); - } + 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); + return (inv_timeScale); - } /* calc_t */ + } /* 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 + /** + * printResJac(): prints out the residual and Jacobian. + * */ #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); + 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(); - printf("\n IntefaceKinetics Object # %d\n\n", isp); + 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); + } - 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"); + 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"); - } - } + 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"); - } - } + 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) { + 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"); - } + 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) { + else if (ioflag > 1) { - updateMFSolnSP(XMolSolnSP); - printf("\n\t Weighted norm of update = %10.4e\n", update_norm); + 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 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"); - } + printf("\t--------------------------------------------------------" + "------------------------------\n"); + } #endif - } /* printIteration */ + } /* 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) { + 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"); - } + 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) { + 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("\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(" 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); + } + 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); + updateMFKinSpecies(XMolKinSpecies, isp); - printf("\n IntefaceKinetics Object # %d\n\n", 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"); + 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"); - } + 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"); - } + 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/spectra/LineBroadener.cpp b/Cantera/src/spectra/LineBroadener.cpp index 4a0e0511d..bfa952de1 100644 --- a/Cantera/src/spectra/LineBroadener.cpp +++ b/Cantera/src/spectra/LineBroadener.cpp @@ -3,25 +3,32 @@ using namespace std; -namespace Cantera { +namespace CanteraSpectra { - doublereal to_wavenumbers(doublereal freq) { - return freq/(100.0*lightSpeed); - } - - doublereal from_wavenumbers(doublereal omega) { - return omega*(100.0*lightSpeed); - } - - Lorentzian::Lorentzian(doublereal FWHM) { - m_hwhm = 0.5*FWHM; + Lorentzian::Lorentzian(doublereal gamma) { + m_hwhm = gamma; m_hwhm2 = m_hwhm*m_hwhm; } + /** + * The Lorentzian profile for collision-broadened lines. + * + *\f[ + * \frac{1}{\pi} \frac{\gamma}{ (\Delta\nu)^2 + \gamma^2} + *\f] + * Units: 1/wavenumber (or cm). + */ doublereal Lorentzian::profile(doublereal deltaFreq) { return (1.0/Pi) *m_hwhm/(deltaFreq*deltaFreq + m_hwhm2); } - + + /** + * + * The cumulative profile, given by + * \f[ + * \frac{1}{\pi} \tan^{-1}\left(\frac{\Delta\nu}{gamma}\right) + 0.5 + * \f] + */ doublereal Lorentzian::cumulative(doublereal deltaFreq) { return (1.0/Pi) * atan(deltaFreq/m_hwhm) + 0.5; } @@ -30,13 +37,15 @@ namespace Cantera { return 2.0*m_hwhm; } - Gaussian::Gaussian(doublereal FWHM) { - m_width = FWHM; - m_sigma = 0.5*FWHM / sqrt(2.0 * log(2.0)); + Gaussian::Gaussian(doublereal sigma) { + m_sigma = sigma; m_sigma2 = m_sigma*m_sigma; } doublereal Gaussian::profile(doublereal deltaFreq) { + //cout << "entered Gaussian::profile" << endl; + //cout << "deltaFreq = " << deltaFreq << endl; + //cout << "m_sigma = " << m_sigma << endl; return 1.0/(m_sigma*SqrtTwo*SqrtPi) * exp(-deltaFreq*deltaFreq/(2.0*m_sigma2)); } @@ -46,7 +55,83 @@ namespace Cantera { } doublereal Gaussian::width() { - return m_width; + return 2.0*m_sigma*sqrt(log(4.0)); + } + + + + /** + * @param sigma The standard deviation of the Gaussian + * @param gamma The half-width of the Lorentzian. + */ + Voigt::Voigt(doublereal sigma, doublereal gamma) { + m_sigma = sigma; + m_sigma2 = m_sigma*m_sigma; + m_gamma_lor = gamma; + m_sigsqrt2 = SqrtTwo*m_sigma; + m_gamma = gamma/m_sigsqrt2; + m_eps = 1.0e-9; + } + + /** + * This method evaluates the function + * \f[ + * F(x, y) = \frac{y}{\pi}\int_{-\infty}^{+\infty} \frac{e^{-z^2}} + * {(x - z)^2 + y^2} dz + * \f] + */ + doublereal Voigt::F(doublereal x) { + + if (x < 0.0) x = -x; + double y = m_gamma; + + double c3 = log(Pi*m_eps/2.0); + double tau = sqrt(-log(y) - c3); + double b = (tau + x)/y; + double t = b*y; + double f1, f2, f3; + double c0 = 2.0/(Pi*m_eps); + const double c1 = 1.0/SqrtTwo; + const double c2 = 2.0/SqrtPi; + + if (y > c0/m_eps) { + throw CanteraError("Voigt::F", + "condition that y < c0/epsion violated"); + } + while (1 > 0) { + f1 = c2*y*exp(-Pi*Pi/(t*t)); + f2 = fabs(y*y - Pi*Pi/(t*t)); + f3 = 1.0 - pow(exp(-Pi*Pi/(t*t)),2); + t *= c1; + // cout << "t = " << t << endl; + if ((f1/(f2*f3)) < 0.5*m_eps) break; + } + double h = t/y; + int N = int(0.5*b/h); + double S = 0.0; + double u = h/2; + for (int i = 0; i < N; i++) { + S += (1.0 + exp(-4.0*x*y*u))*exp(-pow(y*u-x,2))/(u*u+1.0); + u += h; + } + double Q = h*S/Pi; + double C = 0.0; + if (y*y < Pi/h) { + C = 2.0*exp(y*y - x*x)*cos(2*x*y)/(1.0 + exp(2*Pi/h)); + } + else { + return 0.0; + } + return Q + C; + } + + /** + * Voigt profile. + * + */ + doublereal Voigt::profile(doublereal deltaFreq) { + const double ff = m_gamma_lor*m_gamma_lor/(2.0*Pi*m_sigma*m_sigma); + return ff*F(deltaFreq/m_sigsqrt2); } } diff --git a/Cantera/src/spectra/LineBroadener.h b/Cantera/src/spectra/LineBroadener.h index b8420586f..501512449 100644 --- a/Cantera/src/spectra/LineBroadener.h +++ b/Cantera/src/spectra/LineBroadener.h @@ -1,23 +1,28 @@ /** * @file LinerBoadener.h * Header file for class LineBroadener + * @ingroup spectroscopy */ #include "ct_defs.h" #include "ctexceptions.h" -namespace Cantera { +using namespace Cantera; +namespace CanteraSpectra { /** * Base class for classes implementing line shapes of * various types. + * @ingroup spectroscopy */ class LineBroadener { public: + /// Default constructor LineBroadener() {} + /// Destructor virtual ~LineBroadener() {} /** @@ -26,7 +31,8 @@ namespace Cantera { * P(\Delta\nu) *\f] * as a function of distance from line - * center. This function must have total area = 1.0. + * center \f$ \Delta\nu \f$. + * This function must have total area = 1.0. * Note that this method must be overloaded in each * derived class. If the base class method is called, * an exception will be thrown. @@ -36,6 +42,10 @@ namespace Cantera { "base class method called!"); } + doublereal operator()(doublereal deltaFreq) { + return profile(deltaFreq); + } + /** * The cumulative profile, defined as * \f[ @@ -47,6 +57,7 @@ namespace Cantera { "base class method called!"); } + virtual doublereal width() { return 0.0; } }; /** @@ -55,7 +66,7 @@ namespace Cantera { * \f[ * L(\Delta\nu) = \frac{1}{\pi}\frac{\gamma}{\Delta\nu^2 + \gamma^2} * \f] - * where \f$ \gamma = {\mbox{FWHM}/2 \f$. + * where \f$ \gamma = {\mbox{FWHM}/2} \f$. */ class Lorentzian : public LineBroadener { public: @@ -80,15 +91,50 @@ namespace Cantera { * Constructor. * @param FWHM Full width at half-maximum. */ - Gaussian(doublereal FWHM); + Gaussian(doublereal sigma); virtual doublereal profile(doublereal deltaFreq); virtual doublereal cumulative(doublereal deltaFreq); virtual doublereal width(); + doublereal standardDev() { + return m_sigma; + } + protected: doublereal m_sigma; doublereal m_sigma2; doublereal m_width; }; + + /** + * A Voigt profile is the convolution of a Lorentzian and a + * Gaussian profile. This profile results when Doppler + * broadening and collisional broadening both are important. + */ + class Voigt : public LineBroadener { + public: + + /** + * Constructor. + */ + Voigt(doublereal sigma, doublereal gamma); + virtual doublereal profile(doublereal deltaFreq); + //virtual doublereal cumulative(doublereal deltaFreq) + //virtual doublereal width() + + protected: + + doublereal F(doublereal x); + + doublereal m_sigma; + doublereal m_gamma_lor; + doublereal m_sigma2; + doublereal m_width; + doublereal m_gamma; + doublereal m_sigsqrt2; + doublereal m_a; + doublereal m_eps; + }; + } diff --git a/Cantera/src/spectra/rotor.cpp b/Cantera/src/spectra/rotor.cpp index 1aaea1644..7c7385f2e 100644 --- a/Cantera/src/spectra/rotor.cpp +++ b/Cantera/src/spectra/rotor.cpp @@ -10,15 +10,17 @@ #include "ct_defs.h" #include "rotor.h" -namespace Cantera { +namespace CanteraSpectra { /** - * @param Bv Rotational constant, wavenumbers + * Constructor. + * + * @param Bv Rotational constant, wavenumbers. + * @dipoleMoment permanent dipole moment. * @param Dv Coefficient describing centrifugal - * effects on the bond length. For a rigid rotor, Bv = 0. + * effects on the bond length. For a rigid rotor, Bv = 0. * @param Hv Coefficient describing higher-order vibration-rotation * interactions. For a rigid rotor, Hv = 0. - * @dipoleMoment permanent dipole moment. */ Rotor::Rotor(doublereal Bv, doublereal dipoleMoment, doublereal Dv, doublereal Hv ) : m_Bv(Bv), @@ -81,18 +83,23 @@ namespace Cantera { return degeneracy(J)*exp(-wnum_to_J(energy_w(J))/(Boltzmann*T)); } - /** The difference in the energies of an upper and a lower state. - * + /** + * The frequency at which radiation is absorbed by a transition + * from the lower to the upper state in wavenumber units. */ doublereal Rotor::frequency(int J_lower, int J_upper) { return (energy_w(J_upper) - energy_w(J_lower)); } + /** + * The spectral intensity of a rotational transition. + */ doublereal Rotor::intensity(int J_lower, int J_upper, doublereal T) { int dJ = J_upper - J_lower; if (dJ > 1 || dJ < -1) return 0; return relPopulation(J_lower, T); - } + } + } diff --git a/Cantera/src/spectra/rotor.h b/Cantera/src/spectra/rotor.h index 5f27bdb35..70e2ba6ab 100644 --- a/Cantera/src/spectra/rotor.h +++ b/Cantera/src/spectra/rotor.h @@ -1,45 +1,45 @@ +#ifndef CT_ROTOR +#define CT_ROTOR + /** * @file rotor.h * Header file for class Rotor. */ /** - * @defgroup spectra Spectroscopic Models + * @defgroup spectroscopy Spectroscopic Models * * These classes are used to simulate the absorption and emission spectra of * molecules. */ -/* - * @ingroup thermoprops - */ - -#ifndef CT_ROTOR -#define CT_ROTOR - #include "ct_defs.h" +using namespace Cantera; -namespace Cantera { +/** + * Namespace for spectroscopic functions and classes. + */ +namespace CanteraSpectra { /** * Class Rotor represents a non-rigid quantum-mechanical rotor. - * @ingroup spectra + * @ingroup spectroscopy */ class Rotor { public: + + /// Default Constructor. Rotor() {} + + /// Destructor. virtual ~Rotor() {} - /* - */ + /// Full Constructor. Rotor(doublereal Bv, doublereal dipoleMoment = 0.0, doublereal Dv = 0.0, doublereal Hv = 0.0); - - /// energy in wavenumbers doublereal energy_w(int J); - /// degeneracy int degeneracy(int J); doublereal partitionFunction(doublereal T, int cutoff=-1); @@ -66,6 +66,7 @@ namespace Cantera { return freq/(100.0*lightSpeed); } + /** Convert from wavenumbers to Joules. */ inline doublereal wnum_to_J(doublereal w) { return Planck * w * 100.0 * lightSpeed; } diff --git a/Cantera/src/thermo/AdsorbateThermo.h b/Cantera/src/thermo/AdsorbateThermo.h index b09a89d00..1d8230ffd 100644 --- a/Cantera/src/thermo/AdsorbateThermo.h +++ b/Cantera/src/thermo/AdsorbateThermo.h @@ -1,5 +1,5 @@ /** - * @file HarmonicOscThermo.h + * @file AdsorbateThermo.h * * Header for a single-species standard * state object derived from \link Cantera::SpeciesThermoInterpType @@ -21,12 +21,11 @@ #include "SpeciesThermoInterpType.h" -#include -using namespace std; - namespace Cantera { /** + * An adsorbed surface species. + * * This class is designed specifically for use by the class * GeneralSpeciesThermo. It implements a model for the * thermodynamic properties of a molecule that can be modeled as a @@ -86,9 +85,9 @@ namespace Cantera { return (SpeciesThermoInterpType *) np; } - virtual void install(string name, int index, int type, - const doublereal* c, - doublereal minTemp, doublereal maxTemp, doublereal refPressure) { + virtual void install(std::string name, int index, int type, + const doublereal* c, doublereal minTemp, doublereal maxTemp, + doublereal refPressure) { m_be = c[1]; m_nFreqs = int(c[0]); for (int n = 0; n < m_nFreqs; n++) { @@ -120,61 +119,61 @@ namespace Cantera { virtual int speciesIndex() const { return m_index; } - //! Compute the reference-state property of one species - /*! - * Given temperature T in K, this method updates the values of - * the non-dimensional heat capacity at constant pressure, - * enthalpy, and entropy, at the reference pressure, Pref - * of one of the species. The species index is used - * to reference into the cp_R, h_RT, and s_R arrays. - * - * @param temp Temperature (Kelvin) - * @param cp_R Vector of Dimensionless heat capacities. - * (length m_kk). - * @param h_RT Vector of Dimensionless enthalpies. - * (length m_kk). - * @param s_R Vector of Dimensionless entropies. - * (length m_kk). - */ - void updatePropertiesTemp(const doublereal temp, - doublereal* cp_R, - doublereal* h_RT, - doublereal* s_R) const { - h_RT[m_index] = _energy_RT(temp); - cp_R[m_index] = (temp*h_RT[m_index] - - (temp-0.01)*_energy_RT(temp-0.01))/0.01; - s_R[m_index] = h_RT[m_index] - _free_energy_RT(temp); - } - - //!This utility function reports back the type of - //! parameterization and all of the parameters for the - //! species, index. - /*! - * All parameters are output variables - * - * @param n Species index - * @param type Integer type of the standard type - * @param tlow output - Minimum temperature - * @param thigh output - Maximum temperature - * @param pref output - reference pressure (Pa). - * @param coeffs Vector of coefficients used to set the - * parameters for the standard state. - */ - void reportParameters(int &n, int &type, - doublereal &tlow, doublereal &thigh, - doublereal &pref, - doublereal* const coeffs) const { - n = m_index; - type = ADSORBATE; - tlow = m_lowT; - thigh = m_highT; - pref = m_Pref; - coeffs[0] = m_nFreqs; - coeffs[1] = m_be; - for (int i = 2; i < m_nFreqs+2; i++) { - coeffs[i] = m_freq[i-2]; - } - } + //! Compute the reference-state property of one species + /*! + * Given temperature T in K, this method updates the values of + * the non-dimensional heat capacity at constant pressure, + * enthalpy, and entropy, at the reference pressure, Pref + * of one of the species. The species index is used + * to reference into the cp_R, h_RT, and s_R arrays. + * + * @param temp Temperature (Kelvin) + * @param cp_R Vector of Dimensionless heat capacities. + * (length m_kk). + * @param h_RT Vector of Dimensionless enthalpies. + * (length m_kk). + * @param s_R Vector of Dimensionless entropies. + * (length m_kk). + */ + void updatePropertiesTemp(const doublereal temp, + doublereal* cp_R, + doublereal* h_RT, + doublereal* s_R) const { + h_RT[m_index] = _energy_RT(temp); + cp_R[m_index] = (temp*h_RT[m_index] + - (temp-0.01)*_energy_RT(temp-0.01))/0.01; + s_R[m_index] = h_RT[m_index] - _free_energy_RT(temp); + } + + //! This utility function reports back the type of + /*! parameterization and all of the parameters for the + * species, index. + * + * All parameters are output variables + * + * @param n Species index + * @param type Integer type of the standard type + * @param tlow output - Minimum temperature + * @param thigh output - Maximum temperature + * @param pref output - reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the + * parameters for the standard state. + */ + void reportParameters(int &n, int &type, + doublereal &tlow, doublereal &thigh, + doublereal &pref, + doublereal* const coeffs) const { + n = m_index; + type = ADSORBATE; + tlow = m_lowT; + thigh = m_highT; + pref = m_Pref; + coeffs[0] = m_nFreqs; + coeffs[1] = m_be; + for (int i = 2; i < m_nFreqs+2; i++) { + coeffs[i] = m_freq[i-2]; + } + } protected: //! lowest valid temperature @@ -193,36 +192,36 @@ namespace Cantera { doublereal m_be; - doublereal _energy_RT(double T) const { - doublereal x, hnu_kt, hnu, sum = 0.0; - doublereal kt = T*Boltzmann; - int i; - for (i = 0; i < m_nFreqs; i++) { - hnu = Planck * m_freq[i]; - hnu_kt = hnu/kt; - x = exp(-hnu_kt); - sum += hnu_kt * x/(1.0 - x); + doublereal _energy_RT(double T) const { + doublereal x, hnu_kt, hnu, sum = 0.0; + doublereal kt = T*Boltzmann; + int i; + for (i = 0; i < m_nFreqs; i++) { + hnu = Planck * m_freq[i]; + hnu_kt = hnu/kt; + x = exp(-hnu_kt); + sum += hnu_kt * x/(1.0 - x); + } + return sum + m_be/(GasConstant*T); } - return sum + m_be/(GasConstant*T); - } - doublereal _free_energy_RT(double T) const { - doublereal x, hnu_kt, sum = 0.0; - doublereal kt = T*Boltzmann; - int i; - for (i = 0; i < m_nFreqs; i++) { - hnu_kt = Planck * m_freq[i] / kt; - x = exp(-hnu_kt); - sum += log(1.0 - x); + doublereal _free_energy_RT(double T) const { + doublereal x, hnu_kt, sum = 0.0; + doublereal kt = T*Boltzmann; + int i; + for (i = 0; i < m_nFreqs; i++) { + hnu_kt = Planck * m_freq[i] / kt; + x = exp(-hnu_kt); + sum += log(1.0 - x); + } + return sum + m_be/(GasConstant*T); } - return sum + m_be/(GasConstant*T); - } - doublereal _entropy_R(double T) const { - return _energy_RT(T) - _free_energy_RT(T); - } + doublereal _entropy_R(double T) const { + return _energy_RT(T) - _free_energy_RT(T); + } - }; + }; } #endif diff --git a/Cantera/src/thermo/Phase.h b/Cantera/src/thermo/Phase.h index 5935f1773..6b05f1ee4 100755 --- a/Cantera/src/thermo/Phase.h +++ b/Cantera/src/thermo/Phase.h @@ -1,9 +1,12 @@ /** * @file Phase.h - * Header file for class, Phase, which contains functions for setting the - * state of a phase, and for referencing species by name, and also contains text for the module phases - * (see \ref phases and class \link Cantera::Phase Phase\endlink). + * + * Header file for class, Phase, which contains functions for + * setting the state of a phase, and for referencing species by + * name, and also contains text for the module phases (see \ref + * phases and class \link Cantera::Phase Phase\endlink). */ + /* * $Author$ * $Revision$ @@ -25,7 +28,7 @@ namespace Cantera { /** - * @defgroup phases Phases of Matter + * @defgroup phases Models of Phases of Matter * * These classes are used to represent the composition and state of a * single phase of matter. @@ -44,8 +47,9 @@ namespace Cantera { * the phases may then be described using stoichiometry base on the * same Elements class object. * - * The member functions of class %Elements return information about the elements described - * in a particular instantiation of the class. + * The member functions of class %Elements return information about + * the elements described in a particular instantiation of the + * class. * * Class %Constituents is designed to provide information * about the elements and species in a phase - names, index @@ -64,9 +68,9 @@ namespace Cantera { * %Constituents also contains utilities retrieving the index of * a species in the phase given its name, Constituents::speciesIndex(). * - * Class State manages the independent variables of temperature, mass density, - * and species mass/mole fraction that define the thermodynamic - * state. + * Class State manages the independent variables of temperature, + * mass density, and species mass/mole fraction that define the + * thermodynamic state. * * Class %State stores just enough information about a * multicomponent solution to specify its intensive thermodynamic @@ -79,18 +83,20 @@ namespace Cantera { * Class %State is not usually used directly in application * programs. Its primary use is as a base class for class * Phase. Class %State has no virtual methods, and none of its - * methods are meant to be overloaded. However, this is one exception. - * If the phase is incompressible, then the density must be replaced - * by the pressure as the independent variable. In this case, functions - * such as State::setMassFractions() within the class %State must actually now - * calculate the density (at constant T and P) instead of leaving - * it alone as befits an independent variable. Therefore, these types - * of functions are virtual functions and need to be overloaded - * for incompressible phases. Note, for nearly incompressible phases - * (or phases which utilize standard states based on a T and P) this - * change in independent variables may be advantageous as well, - * and these functions in %State need to overload as well so that the - * storred density within State doesn't become out of date. + * methods are meant to be overloaded. However, this is one + * exception. If the phase is incompressible, then the density + * must be replaced by the pressure as the independent variable. In + * this case, functions such as State::setMassFractions() within + * the class %State must actually now calculate the density (at + * constant T and P) instead of leaving it alone as + * befits an independent variable. Therefore, these types of + * functions are virtual functions and need to be overloaded for + * incompressible phases. Note, for nearly incompressible phases + * (or phases which utilize standard states based on a T and + * P) this change in independent variables may be + * advantageous as well, and these functions in %State need to + * overload as well so that the storred density within State + * doesn't become out of date. * * Class Phase derives from both clases * Constituents and State. In addition to the methods of those two diff --git a/Cantera/src/thermo/PureFluidPhase.h b/Cantera/src/thermo/PureFluidPhase.h index 8512a9022..2e2c626c1 100644 --- a/Cantera/src/thermo/PureFluidPhase.h +++ b/Cantera/src/thermo/PureFluidPhase.h @@ -1,12 +1,11 @@ /** * @file PureFluidPhase.h - * Header for a ThermoPhase object for a pure fluid phase consisting of - * gas, liquid, mixed-gas-liquid - * and supercrit fluid (see \ref thermoprops - * and class \link Cantera::PureFluidPhase PureFluidPhase\endlink). - * * - * This object is only available if the WITH_PURE_FLUIDS optional compile + * Header for a ThermoPhase class for a pure fluid phase consisting of + * gas, liquid, mixed-gas-liquid and supercrit fluid (see \ref thermoprops + * and class \link Cantera::PureFluidPhase PureFluidPhase\endlink). + * + * This class is only available if the WITH_PURE_FLUIDS optional compile * capability has been turned on in Cantera's makefile system. * It inherits from ThermoPhase, but is built on top of the tpx package. */ @@ -40,11 +39,12 @@ namespace tpx { namespace Cantera { - //! This phase object consists of a single component that can be a gas, a liquid, - //! a mixed gas-liquid fluid, or a fluid beyond its critical point + //! This phase object consists of a single component that can be a + //! gas, a liquid, a mixed gas-liquid fluid, or a fluid beyond its + //! critical point /*! - * The object inherits from ThermoPhase. However, its build on top of the - * tpx package. + * The object inherits from ThermoPhase. However, its build on top + * of the tpx package. * * *

Specification of Species Standard State Properties

diff --git a/Cantera/src/thermo/SpeciesThermo.h b/Cantera/src/thermo/SpeciesThermo.h index 0989a5190..9c8a5f611 100755 --- a/Cantera/src/thermo/SpeciesThermo.h +++ b/Cantera/src/thermo/SpeciesThermo.h @@ -123,7 +123,7 @@ namespace Cantera { * at a set number of temperatures. Between each temperature * the heat capacity is treated as a constant. * . - * . + * @ingroup phases . */ //@{ diff --git a/Cantera/src/thermo/ThermoPhase.h b/Cantera/src/thermo/ThermoPhase.h index bdfb17e55..9efa7d49c 100755 --- a/Cantera/src/thermo/ThermoPhase.h +++ b/Cantera/src/thermo/ThermoPhase.h @@ -1,8 +1,10 @@ /** * @file ThermoPhase.h + * * Header file for class ThermoPhase, the base class for phases with * thermodynamic properties, and the text for the Module thermoprops - * (see \ref thermoprops and class \link Cantera::ThermoPhase ThermoPhase\endlink). + * (see \ref thermoprops and class \link Cantera::ThermoPhase + * ThermoPhase\endlink). */ /* @@ -174,6 +176,7 @@ namespace Cantera { * read ThermoPhases from XML files. * @see newPhase(XML_Node &phase) How to call the Factory routine to create * and initialize %ThermoPhase objects. + * @ingroup phases */ diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h index c3252d698..7c3ed7f25 100755 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -12,6 +12,14 @@ // Copyright 2001-2003 California Institute of Technology +/** + * @defgroup tranprops Transport Properties + * + * @ingroup phases + * + * These classes provide transport properties. + */ + #ifndef CT_TRANSPORTBASE_H #define CT_TRANSPORTBASE_H