C++ BlackoilState::gasoilratio方法代码示例

        bool equals(const BlackoilState& other, double epsilon = 1e-8) const {
            bool equal = (numPhases() == other.numPhases());

            for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++ phaseIdx) {
                equal = equal && (usedPhases_.phase_used[phaseIdx] == other.usedPhases_.phase_used[phaseIdx]);
                if (usedPhases_.phase_used[phaseIdx])
                    equal = equal && (usedPhases_.phase_pos[phaseIdx] == other.usedPhases_.phase_pos[phaseIdx]);
            }

            equal = equal && (vectorApproxEqual( pressure() , other.pressure() , epsilon));
            equal = equal && (vectorApproxEqual( facepressure() , other.facepressure() , epsilon));
            equal = equal && (vectorApproxEqual( faceflux() , other.faceflux() , epsilon));
            equal = equal && (vectorApproxEqual( surfacevol() , other.surfacevol() , epsilon));
            equal = equal && (vectorApproxEqual( saturation() , other.saturation() , epsilon));
            equal = equal && (vectorApproxEqual( gasoilratio() , other.gasoilratio() , epsilon));

            return equal;
        }

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

        if (pu.phase_used[BlackoilPhases::Vapour]) {
            const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
            phasePressure[gpos][cellIdx] =
                initialState.pressure()[cellIdx]
                + (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]);
        }
    }

    // calculate the densities of the active phases for each cell
    if (pu.phase_used[BlackoilPhases::Aqua]) {
        const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
        const auto& pvtw = props.waterPvt();
        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
            int pvtRegionIdx = pvtRegion[cellIdx];

            double T = initialState.temperature()[cellIdx];
            double p = phasePressure[wpos][cellIdx];
            double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);

            rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
        }
    }

    if (pu.phase_used[BlackoilPhases::Liquid]) {
        const int opos = pu.phase_pos[BlackoilPhases::Liquid];
        const auto& pvto = props.oilPvt();
        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
            int pvtRegionIdx = pvtRegion[cellIdx];

            double T = initialState.temperature()[cellIdx];
            double p = phasePressure[opos][cellIdx];
            double Rs = initialState.gasoilratio()[cellIdx];
            double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p);

            double b;
            if (Rs >= RsSat) {
                b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
            }
            else {
                b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
            }

            rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
            if (pu.phase_used[BlackoilPhases::Vapour]) {
                int gpos = pu.phase_pos[BlackoilPhases::Vapour];
                rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
            }
        }
    }

    if (pu.phase_used[BlackoilPhases::Vapour]) {
        const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
        const auto& pvtg = props.gasPvt();
        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
            int pvtRegionIdx = pvtRegion[cellIdx];

            double T = initialState.temperature()[cellIdx];
            double p = phasePressure[gpos][cellIdx];
            double Rv = initialState.rv()[cellIdx];
            double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p);

            double b;
            if (Rv >= RvSat) {
                b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);

// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
    using namespace Opm;

    std::cout << "\n================    Test program for fully implicit three-phase black-oil flow     ===============\n\n";
    parameter::ParameterGroup param(argc, argv, false);
    std::cout << "---------------    Reading parameters     ---------------" << std::endl;

    // If we have a "deck_filename", grid and props will be read from that.
    bool use_deck = param.has("deck_filename");
    if (!use_deck) {
        OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
                  "Specify the deck with deck_filename=deckname.data (for example).");
    }
    boost::scoped_ptr<EclipseGridParser> deck;
    boost::scoped_ptr<GridManager> grid;
    boost::scoped_ptr<BlackoilPropertiesInterface> props;
    boost::scoped_ptr<BlackoilPropsAdInterface> new_props;
    boost::scoped_ptr<RockCompressibility> rock_comp;
    BlackoilState state;
    // bool check_well_controls = false;
    // int max_well_control_iterations = 0;
    double gravity[3] = { 0.0 };
    std::string deck_filename = param.get<std::string>("deck_filename");
    deck.reset(new EclipseGridParser(deck_filename));
    // Grid init
    grid.reset(new GridManager(*deck));

    // use the capitalized part of the deck's filename between the
    // last '/' and the last '.' character as base name.
    std::string baseName = deck_filename;
    auto charPos = baseName.rfind('/');
    if (charPos != std::string::npos)
        baseName = baseName.substr(charPos + 1);
    charPos = baseName.rfind('.');
    if (charPos != std::string::npos)
        baseName = baseName.substr(0, charPos);
    baseName = boost::to_upper_copy(baseName);

    Opm::EclipseWriter outputWriter(param, share_obj(*deck), share_obj(*grid->c_grid()));
    // Rock and fluid init
    props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param));
    new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid()));
    // check_well_controls = param.getDefault("check_well_controls", false);
    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
    // Rock compressibility.
    rock_comp.reset(new RockCompressibility(*deck));
    // Gravity.
    gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
    // Init state variables (saturation and pressure).
    if (param.has("init_saturation")) {
        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
        initBlackoilSurfvol(*grid->c_grid(), *props, state);
        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
        const PhaseUsage pu = props->phaseUsage();
        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
            const int np = props->numPhases();
            const int nc = grid->c_grid()->number_of_cells;
            for (int c = 0; c < nc; ++c) {
                state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]]
                    / state.surfacevol()[c*np + pu.phase_pos[Oil]];
            }
        }
    } else {
        initBlackoilStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
    }

    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
    const double *grav = use_gravity ? &gravity[0] : 0;

    // Linear solver.
    LinearSolverFactory linsolver(param);

    // Write parameters used for later reference.
    bool output = param.getDefault("output", true);
    std::ofstream epoch_os;
    std::string output_dir;
    if (output) {
        output_dir =
            param.getDefault("output_dir", std::string("output"));
        boost::filesystem::path fpath(output_dir);
        try {
            create_directories(fpath);
        }
        catch (...) {
            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
        }
        std::string filename = output_dir + "/epoch_timing.param";
        epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
        // open file to clean it. The file is appended to in SimulatorTwophase
        filename = output_dir + "/step_timing.param";
        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
        step_os.close();
        param.writeParam(output_dir + "/simulation.param");
    }


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