batchReactor_cvode.cc

Go to the documentation of this file.

#include "batchReactor_cvode.h"
 
 
int rhsf_cvode(sunrealtype t, N_Vector varsCV, N_Vector dvarsdtCV, void *user_data);
 
 
batchReactor_cvode::batchReactor_cvode(std::shared_ptr<Cantera::Solution> cantSol) {
 
    gas = cantSol->thermo(); 
    kin = cantSol->kinetics(); 
 
    nvar = gas->nSpecies();
 
    // Set absolute and relative tolerances
    std::vector<double> atol(nvar, 1E-10);
    double rtol = 1E-4;
 
    // Initialize the integrator with CVODE
    integrator = std::make_shared<integrator_cvode>(rhsf_cvode, this, nvar, rtol, atol);
}
 
 
void batchReactor_cvode::react(double &h, std::vector<double> &y, const double tRun) {
 
    // Store fixed enthalpy and pressure
    h_fixed = h;
    P_fixed = gas->pressure();
 
    // Integrate over the given time
    integrator->integrate(y, tRun);
}
 
 
int batchReactor_cvode::rhsf(const double t, const double *vars, double *dvarsdt) {
 
    // Set mass fractions and state
    gas->setMassFractions_NoNorm(vars);
    gas->setState_HP(h_fixed, P_fixed);
 
    // Calculate density and reaction rates
    rho = gas->density();
    temperature = gas->temperature();
 
    std::vector<double> rr(nvar);
    kin->getNetProductionRates(&rr[0]);
    for (size_t k = 0; k < gas->nSpecies(); k++)
        dvarsdt[k] = rr[k] * gas->molecularWeight(k) / rho;
 
    return 0;
}
 
// CVODE interface; CVODE calls this function, which then calls user_data's rhsf
 
 
int rhsf_cvode(sunrealtype t, N_Vector varsCV, N_Vector dvarsdtCV, void *user_data) {
 
    batchReactor_cvode *bRxr = static_cast<batchReactor_cvode *>(user_data);       // Cast user_data pointer to batchReactor_cvode pointer
 
    double *vars    = N_VGetArrayPointer_Serial(varsCV);                                 // Get array pointers for variables and derivatives
 
    double *dvarsdt = N_VGetArrayPointer_Serial(dvarsdtCV);                              // Call rhsf method of batchReactor_cvode instance
    
    int rv = bRxr->rhsf(t, vars, dvarsdt);
 
    return rv;
}