C++ Kinetics类代码示例

 status_t DLL_EXPORT kin_getdestructionrates_(const integer* n, doublereal* ddot) {
     try {
         Kinetics* k = _fkin(n);
         k->getDestructionRates(ddot);
         return 0;
     }
     catch (CanteraError) {handleError(); return -1;}
 }

 status_t DLL_EXPORT kin_getfwdratesofprogress_(const integer* n, doublereal* fwdROP) {
     Kinetics* k = _fkin(n);
     try {
         k->getFwdRatesOfProgress(fwdROP);
         return 0;
     }
     catch (CanteraError) {handleError(); return -1;}
 }

 status_t DLL_EXPORT kin_getnetratesofprogress_(const integer* n, doublereal* netROP) {
     try {
         Kinetics* k = _fkin(n);
         k->getNetRatesOfProgress(netROP);
         return 0;
     }
     catch (CanteraError) {handleError(); return -1;}
 }

 status_t kin_getnetratesofprogress_(const integer* n, doublereal* netROP)
 {
     try {
         Kinetics* k = _fkin(n);
         k->getNetRatesOfProgress(netROP);
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

 status_t kin_getdestructionrates_(const integer* n, doublereal* ddot)
 {
     try {
         Kinetics* k = _fkin(n);
         k->getDestructionRates(ddot);
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

 status_t kin_getnetproductionrates_(const integer* n, doublereal* wdot)
 {
     try {
         Kinetics* k = _fkin(n);
         k->getNetProductionRates(wdot);
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

 status_t kin_getequilibriumconstants_(const integer* n, doublereal* kc)
 {
     try {
         Kinetics* k = _fkin(n);
         k->getEquilibriumConstants(kc);
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

 status_t DLL_EXPORT kin_getreactionstring_(const integer* n, integer* i, char* buf, ftnlen lenbuf) {
     try {
         Kinetics* k = _fkin(n);
         std::string r = k->reactionString(*i-1);
         int lout = min(lenbuf,r.size());
         std::copy(r.c_str(), r.c_str() + lout, buf);
         for (int nn = lout; nn < lenbuf; nn++) buf[nn] = ' ';
         return 0;
     }
     catch (CanteraError) {handleError(); return -1;}
 }

 status_t DLL_EXPORT kin_advancecoverages_(const integer* n, doublereal* tstep) {
     try {
         Kinetics* k = _fkin(n);
         if (k->type() == cInterfaceKinetics) {
             ((InterfaceKinetics*)k)->advanceCoverages(*tstep);
         }
         else {
             throw CanteraError("kin_advanceCoverages",
                 "wrong kinetics manager type");
         }
         return 0;
     }
     catch (CanteraError) {handleError(); return -1;}
 }

 status_t kin_advancecoverages_(const integer* n, doublereal* tstep)
 {
     try {
         Kinetics* k = _fkin(n);
         if (k->type() == cInterfaceKinetics) {
             ((InterfaceKinetics*)k)->advanceCoverages(*tstep);
         } else {
             throw CanteraError("kin_advanceCoverages",
                                "wrong kinetics manager type");
         }
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

 status_t kin_getreactionstring_(const integer* n, integer* i, char* buf, ftnlen lenbuf)
 {
     try {
         Kinetics* k = _fkin(n);
         std::string r = k->reactionString(*i-1);
         int lout = std::min(lenbuf, (int) r.size());
         std::copy(r.c_str(), r.c_str() + lout, buf);
         for (int nn = lout; nn < lenbuf; nn++) {
             buf[nn] = ' ';
         }
     } catch (...) {
         return handleAllExceptions(-1, ERR);
     }
     return 0;
 }

    status_t ctbuildsolutionfromxml(char* src, integer* ixml, char* id,
                                    integer* ith, integer* ikin, ftnlen lensrc, ftnlen lenid)
    {
        try {
            XML_Node* root = 0;
            if (*ixml > 0) {
                root = _xml(ixml);
            }

            thermo_t* t = _fth(ith);
            Kinetics* k = _fkin(ikin);

            XML_Node* x, *r=0;
            if (root) {
                r = &root->root();
            }
            std::string srcS = f2string(src, lensrc);
            std::string idS  = f2string(id, lenid);
            if (srcS != "") {
                x = get_XML_Node(srcS, r);
            } else {
                x = get_XML_Node(idS, r);
            }
            if (!x) {
                return 0;
            }
            importPhase(*x, t);
            k->addPhase(*t);
            k->init();
            installReactionArrays(*x, *k, x->id());
            t->setState_TP(300.0, OneAtm);
            if (r) {
                if (&x->root() != &r->root()) {
                    delete &x->root();
                }
            } else {
                delete &x->root();
            }
        } catch (...) {
            return handleAllExceptions(-1, ERR);
        }
        return 0;
    }

double Reactor::evalSurfaces(double t, double* ydot)
{
    const vector_fp& mw = m_thermo->molecularWeights();

    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    size_t loc = 0; // offset into ydot
    double mdot_surf = 0.0; // net mass flux from surface

    for (size_t i = 0; i < m_wall.size(); i++) {
        Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
        SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
        if (surf && kin) {
            double rs0 = 1.0/surf->siteDensity();
            size_t nk = surf->nSpecies();
            double sum = 0.0;
            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(&m_work[0]);
            size_t ns = kin->surfacePhaseIndex();
            size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

            double wallarea = m_wall[i]->area();
            for (size_t k = 0; k < m_nsp; k++) {
                m_sdot[k] += m_work[k]*wallarea;
                mdot_surf += m_sdot[k] * mw[k];
            }
        }
    }
    return mdot_surf;
}

bool checkElectrochemReaction(const XML_Node& p, Kinetics& kin, const XML_Node& r)
{
    // If other phases are involved in the reaction mechanism, they must be
    // listed in a 'phaseArray' child element. Homogeneous mechanisms do not
    // need to include a phaseArray element.
    vector<string> phase_ids;
    if (p.hasChild("phaseArray")) {
        const XML_Node& pa = p.child("phaseArray");
        getStringArray(pa, phase_ids);
    }
    phase_ids.push_back(p["id"]);

    // Get reaction product and reactant information
    Composition reactants = parseCompString(r.child("reactants").value());
    Composition products = parseCompString(r.child("products").value());


    // If the reaction has undeclared species don't perform electrochemical check
    for (const auto& sp : reactants) {
        if (kin.kineticsSpeciesIndex(sp.first) == npos) {
            return true;
        }
    }

    for (const auto& sp : products) {
        if (kin.kineticsSpeciesIndex(sp.first) == npos) {
            return true;
        }
    }

    // Initialize the electron counter for each phase
    std::vector<double> e_counter(phase_ids.size(), 0.0);

    // Find the amount of electrons in the products for each phase
    for (const auto& sp : products) {
        const ThermoPhase& ph = kin.speciesPhase(sp.first);
        size_t k = ph.speciesIndex(sp.first);
        double stoich = sp.second;
        for (size_t m = 0; m < phase_ids.size(); m++) {
            if (phase_ids[m] == ph.id()) {
                e_counter[m] += stoich * ph.charge(k);
                break;
            }
        }
    }

    // Subtract the amount of electrons in the reactants for each phase
    for (const auto& sp : reactants) {
        const ThermoPhase& ph = kin.speciesPhase(sp.first);
        size_t k = ph.speciesIndex(sp.first);
        double stoich = sp.second;
        for (size_t m = 0; m < phase_ids.size(); m++) {
            if (phase_ids[m] == ph.id()) {
                e_counter[m] -= stoich * ph.charge(k);
                break;
            }
        }
    }

    // If the electrons change phases then the reaction is electrochemical
    bool echemical = false;
    for(size_t m = 0; m < phase_ids.size(); m++) {
        if (fabs(e_counter[m]) > 1e-4) {
            echemical = true;
            break;
        }
    }

    // If the reaction is electrochemical, ensure the reaction is identified as
    // electrochemical. If not already specified beta is assumed to be 0.5
    std::string type = ba::to_lower_copy(r["type"]);
    if (!r.child("rateCoeff").hasChild("electrochem")) {
        if ((type != "butlervolmer_noactivitycoeffs" &&
             type != "butlervolmer" &&
             type != "surfaceaffinity") &&
             echemical) {
            XML_Node& f = r.child("rateCoeff").addChild("electrochem","");
            f.addAttribute("beta",0.5);
        }
    }
    return true;
}

void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y,
                                   doublereal* ydot, doublereal* params)
{
    size_t nk;
    m_thermo->restoreState(m_state);

    Kinetics* kin;
    size_t npar, ploc;
    double mult;

    // process sensitivity parameters
    if (params) {

        npar = m_pnum.size();
        for (size_t n = 0; n < npar; n++) {
            mult = m_kin->multiplier(m_pnum[n]);
            m_kin->setMultiplier(m_pnum[n], mult*params[n]);
        }
        ploc = npar;
        for (size_t m = 0; m < m_nwalls; m++) {
            if (m_nsens_wall[m] > 0) {
                m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
                ploc += m_nsens_wall[m];
            }
        }
    }

    m_Q = 0.0;

    // compute wall terms
    doublereal rs0, sum, wallarea;
    double mcpdTdt = 0.0; // m * c_p * dT/dt
    double dmdt = 0.0; // dm/dt (gas phase)
    double* dYdt = ydot + 2;
    m_thermo->getPartialMolarEnthalpies(&m_hk[0]);

    SurfPhase* surf;
    size_t lr, ns, loc = m_nsp+2, surfloc;
    fill(m_sdot.begin(), m_sdot.end(), 0.0);
    for (size_t i = 0; i < m_nwalls; i++) {
        lr = 1 - 2*m_lr[i];
        m_Q += lr*m_wall[i]->Q(time);
        kin = m_wall[i]->kinetics(m_lr[i]);
        surf = m_wall[i]->surface(m_lr[i]);
        if (surf && kin) {
            rs0 = 1.0/surf->siteDensity();
            nk = surf->nSpecies();
            sum = 0.0;
            surf->setTemperature(m_state[0]);
            m_wall[i]->syncCoverages(m_lr[i]);
            kin->getNetProductionRates(DATA_PTR(m_work));
            ns = kin->surfacePhaseIndex();
            surfloc = kin->kineticsSpeciesIndex(0,ns);
            for (size_t k = 1; k < nk; k++) {
                ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
                sum -= ydot[loc + k];
            }
            ydot[loc] = sum;
            loc += nk;

            wallarea = m_wall[i]->area();
            for (size_t k = 0; k < m_nsp; k++) {
                m_sdot[k] += m_work[k]*wallarea;
            }
        }
    }

    const vector_fp& mw = m_thermo->molecularWeights();
    const doublereal* Y = m_thermo->massFractions();

    if (m_chem) {
        m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot"
    }

    double mdot_surf = 0.0; // net mass flux from surface
    for (size_t k = 0; k < m_nsp; k++) {
        // production in gas phase and from surfaces
        dYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k] / m_mass;
        mdot_surf += m_sdot[k] * mw[k];
    }
    dmdt += mdot_surf;

    // external heat transfer
    mcpdTdt -= m_Q;

    for (size_t n = 0; n < m_nsp; n++) {
        // heat release from gas phase and surface reations
        mcpdTdt -= m_wdot[n] * m_hk[n] * m_vol;
        mcpdTdt -= m_sdot[n] * m_hk[n];
        // dilution by net surface mass flux
        dYdt[n] -= Y[n] * mdot_surf / m_mass;
    }

    // add terms for open system
    if (m_open) {
        // outlets
        for (size_t i = 0; i < m_nOutlets; i++) {
            dmdt -= m_outlet[i]->massFlowRate(time); // mass flow out of system
        }

//.........这里部分代码省略.........