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

 /// @brief Computes injected and produced surface volumes of all phases.
 /// Note 1: assumes that only the first phase is injected.
 /// Note 2: assumes that transport has been done with an
 ///         implicit method, i.e. that the current state
 ///         gives the mobilities used for the preceding timestep.
 /// Note 3: Gives surface volume values, not reservoir volumes
 ///         (as the incompressible version of the function does).
 ///         Also, assumes that transport_src is given in surface volumes
 ///         for injector terms!
 /// @param[in]  props           fluid and rock properties.
 /// @param[in]  state           state variables (pressure, sat, surfvol)
 /// @param[in]  transport_src   if < 0: total resv outflow, if > 0: first phase surfv inflow
 /// @param[in]  dt              timestep used
 /// @param[out] injected        must point to a valid array with P elements,
 ///                             where P = s.size()/src.size().
 /// @param[out] produced        must also point to a valid array with P elements.
 void computeInjectedProduced(const BlackoilPropertiesInterface& props,
                              const BlackoilState& state,
                              const std::vector<double>& transport_src,
                              const double dt,
                              double* injected,
                              double* produced)
 {
     const int num_cells = transport_src.size();
     if (props.numCells() != num_cells) {
         OPM_THROW(std::runtime_error, "Size of transport_src vector does not match number of cells in props.");
     }
     const int np = props.numPhases();
     if (int(state.saturation().size()) != num_cells*np) {
         OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
     }
     const std::vector<double>& press = state.pressure();
     const std::vector<double>& temp = state.temperature();
     const std::vector<double>& s = state.saturation();
     const std::vector<double>& z = state.surfacevol();
     std::fill(injected, injected + np, 0.0);
     std::fill(produced, produced + np, 0.0);
     std::vector<double> visc(np);
     std::vector<double> mob(np);
     std::vector<double> A(np*np);
     std::vector<double> prod_resv_phase(np);
     std::vector<double> prod_surfvol(np);
     for (int c = 0; c < num_cells; ++c) {
         if (transport_src[c] > 0.0) {
             // Inflowing transport source is a surface volume flux
             // for the first phase.
             injected[0] += transport_src[c]*dt;
         } else if (transport_src[c] < 0.0) {
             // Outflowing transport source is a total reservoir
             // volume flux.
             const double flux = -transport_src[c]*dt;
             const double* sat = &s[np*c];
             props.relperm(1, sat, &c, &mob[0], 0);
             props.viscosity(1, &press[c], &temp[c], &z[np*c], &c, &visc[0], 0);
             props.matrix(1, &press[c], &temp[c], &z[np*c], &c, &A[0], 0);
             double totmob = 0.0;
             for (int p = 0; p < np; ++p) {
                 mob[p] /= visc[p];
                 totmob += mob[p];
             }
             std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0);
             for (int p = 0; p < np; ++p) {
                 prod_resv_phase[p] = (mob[p]/totmob)*flux;
                 for (int q = 0; q < np; ++q) {
                     prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p];
                 }
             }
             for (int p = 0; p < np; ++p) {
                 produced[p] += prod_surfvol[p];
             }
         }
     }
 }

    /// @brief Computes saturation from surface volume
    void computeSaturation(const BlackoilPropertiesInterface& props,
                           BlackoilState& state)
    {

        const int np = props.numPhases();
        const int nc = props.numCells();
        std::vector<double> allA(nc*np*np);
        std::vector<int> allcells(nc);
        for (int c = 0; c < nc; ++c) {
            allcells[c] = c;
        }

        //std::vector<double> res_vol(np);
        const std::vector<double>& z = state.surfacevol();

        props.matrix(nc, &state.pressure()[0], &z[0], &allcells[0], &allA[0], 0);

        // Linear solver.
        MAT_SIZE_T n = np;
        MAT_SIZE_T nrhs = 1;
        MAT_SIZE_T lda = np;
        std::vector<MAT_SIZE_T> piv(np);
        MAT_SIZE_T ldb = np;
        MAT_SIZE_T info = 0;


        //double res_vol;
        double tot_sat;
        const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());

        for (int c = 0; c < nc; ++c) {
            double* A = &allA[c*np*np];
            const double* z_loc = &z[c*np];
            double* s = &state.saturation()[c*np];

            for (int p = 0; p < np; ++p){
                s[p] = z_loc[p];
                }

            dgesv_(&n, &nrhs, &A[0], &lda, &piv[0], &s[0], &ldb, &info);

            tot_sat = 0;
            for (int p = 0; p < np; ++p){
                if (s[p] < epsilon) // saturation may be less then zero due to round of errors
                    s[p] = 0;

                tot_sat += s[p];
            }

            for (int p = 0; p < np; ++p){
                s[p]  = s[p]/tot_sat;
            }





        }

    }

 /// Compute per-solve dynamic properties.
 void CompressibleTpfaPolymer::computePerSolveDynamicData(const double /* dt */,
                                                          const BlackoilState& state,
                                                          const WellState& /* well_state */)
 {
     // std::vector<double> cell_relperm__;
     // std::vector<double> cell_eff_relperm_;
     const int nc = grid_.number_of_cells;
     const int np = props_.numPhases();
     cell_relperm_.resize(nc*np);
     cell_eff_viscosity_.resize(nc*np);
     const double* cell_s = &state.saturation()[0];
     props_.relperm(nc, cell_s, &allcells_[0], &cell_relperm_[0], 0);
     computeWellPotentials(state);
     if (rock_comp_props_ && rock_comp_props_->isActive()) {
         computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
     }
 }

    /// Compute per-iteration dynamic properties for cells.
    void CompressibleTpfa::computeCellDynamicData(const double /*dt*/,
                                                  const BlackoilState& state,
                                                  const WellState& /*well_state*/)
    {
        // These are the variables that get computed by this function:
        //
        // std::vector<double> cell_A_;
        // std::vector<double> cell_dA_;
        // std::vector<double> cell_viscosity_;
        // std::vector<double> cell_phasemob_;
        // std::vector<double> cell_voldisc_;
        // std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
        // std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
        const int nc = grid_.number_of_cells;
        const int np = props_.numPhases();
        const double* cell_p = &state.pressure()[0];
        const double* cell_z = &state.surfacevol()[0];
        const double* cell_s = &state.saturation()[0];
        cell_A_.resize(nc*np*np);
        cell_dA_.resize(nc*np*np);
        props_.matrix(nc, cell_p, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]);
        cell_viscosity_.resize(nc*np);
        props_.viscosity(nc, cell_p, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
        cell_phasemob_.resize(nc*np);
        props_.relperm(nc, cell_s, &allcells_[0], &cell_phasemob_[0], 0);
        std::transform(cell_phasemob_.begin(), cell_phasemob_.end(),
                       cell_viscosity_.begin(),
                       cell_phasemob_.begin(),
                       std::divides<double>());
        // Volume discrepancy: we have that
        //     z = Au, voldiscr = sum(u) - 1,
        // but I am not sure it is actually needed.
        // Use zero for now.
        // TODO: Check this!
        cell_voldisc_.clear();
        cell_voldisc_.resize(nc, 0.0);

        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_);
            rock_comp_.resize(nc);
            for (int cell = 0; cell < nc; ++cell) {
                rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]);
            }
        }
    }

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

    /// Compute well potentials.
    void CompressibleTpfa::computeWellPotentials(const BlackoilState& state)
    {
        if (wells_ == NULL) return;

        const int nw = wells_->number_of_wells;
        const int np = props_.numPhases();
        const int nperf = wells_->well_connpos[nw];
        const int dim = grid_.dimensions;
        const double grav = gravity_ ? gravity_[dim - 1] : 0.0;
        wellperf_wdp_.clear();
        wellperf_wdp_.resize(nperf, 0.0);
        if (not (std::abs(grav) > 0.0)) {
            return;
        }

        // Temporary storage for perforation A matrices and densities.
        std::vector<double> A(np*np, 0.0);
        std::vector<double> rho(np, 0.0);

        // Main loop, iterate over all perforations,
        // using the following formula (by phase):
        //    wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
        // where the total density rho(perf) is taken to be
        //    sum_p (rho_p*saturation_p) in the perforation cell.
        for (int w = 0; w < nw; ++w) {
            const double ref_depth = wells_->depth_ref[w];
            for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w + 1]; ++j) {
                const int cell = wells_->well_cells[j];
                const double cell_depth = grid_.cell_centroids[dim * cell + dim - 1];
                props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
                props_.density(1, &A[0], &cell, &rho[0]);
                for (int phase = 0; phase < np; ++phase) {
                    const double s_phase = state.saturation()[np*cell + phase];
                    wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
                }
            }
        }
    }

    SimulatorReport SimulatorCompressibleTwophase::Impl::run(SimulatorTimer& timer,
                                                             BlackoilState& state,
                                                             WellState& well_state)
    {
        std::vector<double> transport_src;

        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        std::vector<double> initial_porevol = porevol;

        // Main simulation loop.
        Opm::time::StopWatch pressure_timer;
        double ptime = 0.0;
        Opm::time::StopWatch transport_timer;
        double ttime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        double init_surfvol[2] = { 0.0 };
        double inplace_surfvol[2] = { 0.0 };
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
            wellreport.push(props_, *wells_,
                            state.pressure(), state.surfacevol(), state.saturation(),
                            0.0, well_state.bhp(), well_state.perfRates());
        }
        std::fstream tstep_os;
        if (output_) {
            std::string filename = output_dir_ + "/step_timing.param";
            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
        }
        for (; !timer.done(); ++timer) {
            // Report timestep and (optionally) write state to disk.
            step_timer.start();
            timer.report(std::cout);
            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                if (output_vtk_) {
                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                }
                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
            }

            SimulatorReport sreport;

            // Solve pressure equation.
            if (check_well_controls_) {
                computeFractionalFlow(props_, allcells_,
                                      state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
                                      fractional_flows);
                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
            }
            bool well_control_passed = !check_well_controls_;
            int well_control_iteration = 0;
            do {
                // Run solver.
                pressure_timer.start();
                std::vector<double> initial_pressure = state.pressure();
                psolver_.solve(timer.currentStepLength(), state, well_state);

                // Renormalize pressure if both fluids and rock are
                // incompressible, and there are no pressure
                // conditions (bcs or wells).  It is deemed sufficient
                // for now to renormalize using geometric volume
                // instead of pore volume.
                if (psolver_.singularPressure()) {
                    // Compute average pressures of previous and last
                    // step, and total volume.
                    double av_prev_press = 0.0;
                    double av_press = 0.0;
                    double tot_vol = 0.0;
                    const int num_cells = grid_.number_of_cells;
                    for (int cell = 0; cell < num_cells; ++cell) {
                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
                        tot_vol       += grid_.cell_volumes[cell];
                    }
                    // Renormalization constant
                    const double ren_const = (av_prev_press - av_press)/tot_vol;
                    for (int cell = 0; cell < num_cells; ++cell) {
                        state.pressure()[cell] += ren_const;
                    }
                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
                    for (int well = 0; well < num_wells; ++well) {
                        well_state.bhp()[well] += ren_const;
                    }
                }
//.........这里部分代码省略.........

    SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
                                                              BlackoilState& state,
                                                              WellState& well_state)
    {
        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        std::vector<double> initial_porevol = porevol;

        // Main simulation loop.
        Opm::time::StopWatch solver_timer;
        double stime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
#if 0
        // These must be changed for three-phase.
        double init_surfvol[2] = { 0.0 };
        double inplace_surfvol[2] = { 0.0 };
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
#endif
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
#if 0
            wellreport.push(props_, *wells_,
                            state.pressure(), state.surfacevol(), state.saturation(),
                            0.0, well_state.bhp(), well_state.perfRates());
#endif
        }
        std::fstream tstep_os;
        if (output_) {
            std::string filename = output_dir_ + "/step_timing.param";
            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
        }
        for (; !timer.done(); ++timer) {
            // Report timestep and (optionally) write state to disk.
            step_timer.start();
            timer.report(std::cout);
            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                if (output_vtk_) {
                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                }
                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);

            }

            SimulatorReport sreport;

            // Solve pressure equation.
            // if (check_well_controls_) {
            //     computeFractionalFlow(props_, allcells_,
            //                           state.pressure(), state.surfacevol(), state.saturation(),
            //                           fractional_flows);
            //     wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
            // }
            bool well_control_passed = !check_well_controls_;
            int well_control_iteration = 0;
            do {
                // Run solver.
                solver_timer.start();
                std::vector<double> initial_pressure = state.pressure();
                solver_.step(timer.currentStepLength(), state, well_state);

                // Stop timer and report.
                solver_timer.stop();
                const double st = solver_timer.secsSinceStart();
                std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;
                stime += st;
                sreport.pressure_time = st;

                // Optionally, check if well controls are satisfied.
                if (check_well_controls_) {
                    Opm::computePhaseFlowRatesPerWell(*wells_,
                                                      well_state.perfRates(),
                                                      fractional_flows,
                                                      well_resflows_phase);
                    std::cout << "Checking well conditions." << std::endl;
                    // For testing we set surface := reservoir
                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
                    ++well_control_iteration;
                    if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
                        OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
                    }
                    if (!well_control_passed) {
                        std::cout << "Well controls not passed, solving again." << std::endl;
                    } else {
                        std::cout << "Well conditions met." << std::endl;
//.........这里部分代码省略.........

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;
    }
//.........这里部分代码省略.........