Method eval

src/zeroD/Reactor.cpp:219–281

}

void Reactor::eval(double time, double* LHS, double* RHS)
{
    double& dmdt = RHS[0];
    double* mdYdt = RHS + 3; // mass * dY/dt

    evalWalls(time);
    m_thermo->restoreState(m_state);
    const vector<double>& mw = m_thermo->molecularWeights();
    const double* Y = m_thermo->massFractions();

    evalSurfaces(LHS + m_nsp + 3, RHS + m_nsp + 3, m_sdot.data());
     // mass added to gas phase from surface reactions
    double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin());
    dmdt = mdot_surf;

    // volume equation
    RHS[1] = m_vdot;

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

    for (size_t k = 0; k < m_nsp; k++) {
        // production in gas phase and from surfaces
        mdYdt[k] = (m_wdot[k] * m_vol + m_sdot[k]) * mw[k];
        // dilution by net surface mass flux
        mdYdt[k] -= Y[k] * mdot_surf;
        LHS[k+3] = m_mass;
    }

    // Energy equation.
    // @f[
    //     \dot U = -P\dot V + A \dot q + \dot m_{in} h_{in} - \dot m_{out} h.
    // @f]
    if (m_energy) {
        RHS[2] = - m_thermo->pressure() * m_vdot + m_Qdot;
    } else {
        RHS[2] = 0.0;
    }

    // add terms for outlets
    for (auto outlet : m_outlet) {
        double mdot = outlet->massFlowRate();
        dmdt -= mdot; // mass flow out of system
        if (m_energy) {
            RHS[2] -= mdot * m_enthalpy;
        }
    }

    // add terms for inlets
    for (auto inlet : m_inlet) {
        double mdot = inlet->massFlowRate();
        dmdt += mdot; // mass flow into system
        for (size_t n = 0; n < m_nsp; n++) {
            double mdot_spec = inlet->outletSpeciesMassFlowRate(n);
            // flow of species into system and dilution by other species
            mdYdt[n] += (mdot_spec - mdot * Y[n]);
        }