twoTanksMacroscopicProblemCxx.C

Go to the documentation of this file.

#include <stdio.h>
#include <iostream>
#include "modenaModel.H"

using namespace std;

int
main(int argc, char *argv[])
{
    const double D = 0.01;

    double p0 = 3e5;
    double p1 = 10000;
    double V0 = 0.1;
    double V1 = 1;
    double T = 300;

    double t = 0.0;
    const double deltat = 1e-3;
    const double tend = 5.5;

    double m0 = p0*V0/287.1/T;
    double m1 = p1*V1/287.1/T;

    double rho0 = m0/V0;
    double rho1 = m1/V1;

    try
    {
        // Instantiate a model
        Modena::modenaModel model("flowRate");

        cout << "inputs:" << endl;
        std::vector<std::string> iNames = model.inputs_names();
        for(int i=0; i<model.inputs_size(); i++)
        {
            cout << iNames[i] << endl;
        }

        cout << "outputs:" << endl;
        std::vector<std::string> oNames = model.outputs_names();
        for(int i=0; i<model.outputs_size(); i++)
        {
            cout << oNames[i] << endl;
        }

        cout << "parameters:" << endl;
        std::vector<std::string> pNames = model.parameters_names();
        for(int i=0; i<model.parameters_size(); i++)
        {
            cout << pNames[i] << endl;
        }

        size_t Dpos = model.inputs_argPos("D");
        size_t rho0Pos = model.inputs_argPos("rho0");
        size_t p0Pos = model.inputs_argPos("p0");
        size_t p1Byp0Pos = model.inputs_argPos("p1Byp0");

        model.argPos_check();

        while(t + deltat < tend + 1e-10)
        {
            t += deltat;

            if(p0 > p1)
            {
                // Set input vector
                model.inputs_set(Dpos, D);
                model.inputs_set(rho0Pos, rho0);
                model.inputs_set(p0Pos, p0);
                model.inputs_set(p1Byp0Pos, p1/p0);

                // Call the model
                model.call();

                // Fetch result
                double mdot = model.outputs_get(0);

                m0 -= mdot*deltat;
                m1 += mdot*deltat;
            }
            else
            {
                // Set input vector
                model.inputs_set(Dpos, D);
                model.inputs_set(rho0Pos, rho1);
                model.inputs_set(p0Pos, p1);
                model.inputs_set(p1Byp0Pos, p0/p1);

                // Call the model
                model.call();

                // Fetch result
                double mdot = model.outputs_get(0);

                m0 += mdot*deltat;
                m1 -= mdot*deltat;
            }

            rho0 = m0/V0;
            rho1 = m1/V1;
            p0 = m0/V0*287.1*T;
            p1 = m1/V1*287.1*T;

            cout << "t = " << t << " rho0 = " << rho0 << " p0 = " << p0 << " p1 = " << p1 << endl;
        }
    }
    catch(const Modena::modenaException& e)
    {
        return e.errorCode();
    }

    return 0;
}