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

 /// \brief Constructor.
 /// \param sendGrid   The grid that the data is attached to when sending.
 /// \param recvGrid   The grid that the data is attached to when receiving.
 /// \param sendState  The state where we will retieve the values to be sent.
 /// \param recvState The state where we will store the received values.
 BlackoilStateDataHandle(const Dune::CpGrid& sendGrid,
                         const Dune::CpGrid& recvGrid,
                         const BlackoilState& sendState,
                         BlackoilState& recvState)
     : sendGrid_(sendGrid), recvGrid_(recvGrid), sendState_(sendState), recvState_(recvState)
 {
     // construction does not resize surfacevol and hydroCarbonState. Do it manually.
     recvState.surfacevol().resize(recvGrid.numCells()*sendState.numPhases(),
                                   std::numeric_limits<double>::max());
     recvState.hydroCarbonState().resize(recvGrid.numCells());
 }

        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;
        }

void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
                  const DeckConstPtr& deck,
                  EclipseStateConstPtr eclipseState,
                  const Grid& grid,
                  const BlackoilState& initialState,
                  const BlackoilPropertiesFromDeck& props,
                  const double gravity)
{

    const PhaseUsage& pu = props.phaseUsage();

    const auto& eqlnum = eclipseState->get3DProperties().getIntGridProperty("EQLNUM");
    const auto& eqlnumData = eqlnum.getData();

    const int numPhases = initialState.numPhases();
    const int numCells = UgGridHelpers::numCells(grid);
    const int numPvtRegions = deck->getKeyword("TABDIMS").getRecord(0).getItem("NTPVT").get< int >(0);

    // retrieve the minimum (residual!?) and the maximum saturations for all cells
    std::vector<double> minSat(numPhases*numCells);
    std::vector<double> maxSat(numPhases*numCells);
    std::vector<int> allCells(numCells);
    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
        allCells[cellIdx] = cellIdx;
    }
    props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());

    // retrieve the surface densities
    std::vector<std::vector<double> > surfaceDensity(numPvtRegions);
    const auto& densityKw = deck->getKeyword("DENSITY");
    for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) {
        surfaceDensity[regionIdx].resize(numPhases);

        if (pu.phase_used[BlackoilPhases::Aqua]) {
            const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
            surfaceDensity[regionIdx][wpos] =
                densityKw.getRecord(regionIdx).getItem("WATER").getSIDouble(0);
        }

        if (pu.phase_used[BlackoilPhases::Liquid]) {
            const int opos = pu.phase_pos[BlackoilPhases::Liquid];
            surfaceDensity[regionIdx][opos] =
                densityKw.getRecord(regionIdx).getItem("OIL").getSIDouble(0);
        }

        if (pu.phase_used[BlackoilPhases::Vapour]) {
            const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
            surfaceDensity[regionIdx][gpos] =
                densityKw.getRecord(regionIdx).getItem("GAS").getSIDouble(0);
        }
    }

    // retrieve the PVT region of each cell. note that we need c++ instead of
    // Fortran indices.
    const int* gc = UgGridHelpers::globalCell(grid);
    std::vector<int> pvtRegion(numCells);
    const auto& cartPvtRegion = eclipseState->get3DProperties().getIntGridProperty("PVTNUM").getData();
    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
        const int cartCellIdx = gc ? gc[cellIdx] : cellIdx;
        pvtRegion[cellIdx] = std::max(0, cartPvtRegion[cartCellIdx] - 1);
    }

    // compute the initial "phase presence" of each cell (required to calculate
    // the inverse formation volume factors
    std::vector<PhasePresence> cond(numCells);
    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
        if (pu.phase_used[BlackoilPhases::Aqua]) {
            const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
            if (sw > 0.0) {
                cond[cellIdx].setFreeWater();
            }
        }

        if (pu.phase_used[BlackoilPhases::Liquid]) {
            const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
            if (so > 0.0) {
                cond[cellIdx].setFreeOil();
            }
        }

        if (pu.phase_used[BlackoilPhases::Vapour]) {
            const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
            if (sg > 0.0) {
                cond[cellIdx].setFreeGas();
            }
        }
    }

    // calculate the initial fluid densities for the gravity correction.
    std::vector<std::vector<double>> rho(numPhases);
    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
        rho[phaseIdx].resize(numCells);
    }

    // compute the capillary pressures of the active phases
    std::vector<double> capPress(numCells*numPhases);
    std::vector<int> cellIdxArray(numCells);
    for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
        cellIdxArray[cellIdx] = cellIdx;
    }
//.........这里部分代码省略.........

inline
void distributeGridAndData( Dune::CpGrid& grid,
                            Opm::DeckConstPtr deck,
                            EclipseStateConstPtr eclipseState,
                            BlackoilState& state,
                            BlackoilPropsAdFromDeck& properties,
                            DerivedGeology& geology,
                            std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>& material_law_manager,
                            std::vector<double>& threshold_pressures,
                            boost::any& parallelInformation,
                            const bool useLocalPerm)
{
    Dune::CpGrid global_grid ( grid );
    global_grid.switchToGlobalView();

    // distribute the grid and switch to the distributed view
    grid.loadBalance(eclipseState, geology.transmissibility().data());
    grid.switchToDistributedView();
    std::vector<int> compressedToCartesianIdx;
    Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
    typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
    auto distributed_material_law_manager = std::make_shared<MaterialLawManager>();
    distributed_material_law_manager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
    // copy the values from the global to the local MaterialLawManager
    // We should actually communicate these to be future proof. But that is
    // really, really cumbersome for the underlying vector<shared_ptr>
    // where the classes pointed to even have more shared_ptr stored in them.
    typedef Dune::CpGrid::ParallelIndexSet IndexSet;
    const IndexSet& local_indices  = grid.getCellIndexSet();
    for ( auto index : local_indices )
    {
        distributed_material_law_manager->materialLawParamsPointerReferenceHack(index.local()) =
            material_law_manager->materialLawParamsPointerReferenceHack(index.global());

        distributed_material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.local()) =
            material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.global());
    }
    BlackoilPropsAdFromDeck distributed_props(properties,
                                              distributed_material_law_manager,
                                              grid.numCells());
    BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
    BlackoilStateDataHandle state_handle(global_grid, grid,
                                         state, distributed_state);
    BlackoilPropsDataHandle props_handle(properties,
                                         distributed_props);
    grid.scatterData(state_handle);
    grid.scatterData(props_handle);
    // Create a distributed Geology. Some values will be updated using communication
    // below
    DerivedGeology distributed_geology(grid,
                                       distributed_props, eclipseState,
                                       useLocalPerm, geology.gravity());
    GeologyDataHandle geo_handle(global_grid, grid,
                                 geology, distributed_geology);
    grid.scatterData(geo_handle);

    std::vector<double> distributed_pressures;

    if( !threshold_pressures.empty() ) // Might be empty if not specified
    {
        if( threshold_pressures.size() !=
            static_cast<std::size_t>(UgGridHelpers::numFaces(global_grid)) )
        {
            OPM_THROW(std::runtime_error, "NNCs not yet supported for parallel runs. "
                      << UgGridHelpers::numFaces(grid) << " faces but " <<
                      threshold_pressures.size()<<" threshold pressure values");
        }
        distributed_pressures.resize(UgGridHelpers::numFaces(grid));
        ThresholdPressureDataHandle press_handle(global_grid, grid,
                                                 threshold_pressures,
                                                 distributed_pressures);
        grid.scatterData(press_handle);
    }

    // copy states
    properties           = distributed_props;
    geology              = distributed_geology;
    state                = distributed_state;
    material_law_manager = distributed_material_law_manager;
    threshold_pressures   = distributed_pressures;
    extractParallelGridInformationToISTL(grid, parallelInformation);
}