Cantera reactor network restart example

reactornet-restart.cpp

// Example demonstrating how to restart a Reactor network using a new
// initial condition.

#include <cantera/zerodim.h>
#include <cantera/IdealGasMix.h>

#include <iostream>

using namespace Cantera;

void runexample()
{
    // use reaction mechanism GRI-Mech 3.0
    IdealGasMix gas("gri30.cti", "gri30");
    size_t nsp = gas.nSpecies();

    // create the combustor, and fill it in initially with N2
    gas.setState_TPX(850.0, OneAtm,
                     "CO:1.0, O2:0.2, CO2:0.5, H2O:0.2, H2:0.1, N2:4");
    Reactor combustor;
    combustor.insert(gas);
    combustor.setInitialVolume(1.0);

    // the simulation only contains one reactor
    ReactorNet sim;
    sim.addReactor(&combustor);

    // A place to save values halfway through the integration
    double T_mid, P_mid, t_mid;
    std::vector<double> Y_mid(nsp);
    bool savedMid = false;

    // Integrate
    double tfinal = 1.0;
    double dt = 0.05;

    std::ofstream f("combustor-restart.csv");

    for (double t=dt; t < tfinal; t += dt) {
        sim.advance(t);

        ThermoPhase& c = combustor.contents();

        if (!savedMid && t >= tfinal/2) {
            savedMid = true;
            T_mid = c.temperature();
            P_mid = c.pressure();
            t_mid = t;
            for (size_t k = 0; k < nsp; k++) {
                Y_mid[k] = c.massFraction(k);
            }
        }

        f << t << ", "
          << combustor.temperature() << ", ";
        for (size_t k = 0; k < nsp; k++) {
            f << c.moleFraction(k);
            if (k != nsp-1) {
                f << ", ";
            }
        }
        f << std::endl;
    }
    f.close();

    // Reset the reactor state to t = tfinal/2
    // If you have a network with more than one reactor,
    // You will need to set the state for each of them.
    gas.setState_TPY(T_mid, P_mid, &Y_mid[0]);
    combustor.insert(gas);

    // Do the integration again, but only from the halfway point.
    sim.initialize(t_mid);

    std::ofstream f2("combustor-restart2.csv");

    for (double t=t_mid + dt; t < tfinal; t += dt) {
        sim.advance(t);

        ThermoPhase& c = combustor.contents();
        f2 << t << ", "
          << combustor.temperature() << ", ";
        for (size_t k = 0; k < nsp; k++) {
            f2 << c.moleFraction(k);
            if (k != nsp-1) {
                f2 << ", ";
            }
        }
        f2 << std::endl;
    }
    f2.close();

    // If you compare the two output files, you will see that they
    // agree to within the integrator's tolerances for the period
    // t_mid < t < tfinal.
}

int main(int argc, char** argv)
{
    try {
        runexample();
        return 0;
    } catch (CanteraError) {
        showErrors(std::cout);
        std::cout << " terminating... " << std::endl;
        appdelete();
        return 1;
    }
}