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

    /// @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], &state.temperature()[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;
            }





        }

    }

 /// @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];
             }
         }
     }
 }

            /**
             * Compute average temperature in all regions.
             *
             * \param[in] state Dynamic reservoir state.
             */
            void
            averageTemperature(const BlackoilState& state)
            {
                T_avg_.setZero();

                const std::vector<double>& T = state.temperature();
                for (std::vector<double>::size_type
                         i = 0, n = T.size(); i < n; ++i)
                {
                    T_avg_(rmap_.region(i)) += T[i];
                }

                T_avg_ /= ncells_;
            }

    /// Compute per-iteration dynamic properties for cells.
    void CompressibleTpfaPolymer::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_eff_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_T = &state.temperature()[0];
        const double* cell_z = &state.surfacevol()[0];
        cell_A_.resize(nc*np*np);
        cell_dA_.resize(nc*np*np);
        props_.matrix(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_A_[0], &cell_dA_[0]);
        cell_viscosity_.resize(nc*np);
        props_.viscosity(nc, cell_p, cell_T, cell_z, &allcells_[0], &cell_viscosity_[0], 0);
        cell_phasemob_.resize(nc*np);
        for (int cell = 0; cell < nc; ++cell) {
            poly_props_.effectiveVisc((*c_)[cell], cell_viscosity_[np*cell + 0], cell_eff_viscosity_[np*cell + 0]);
            poly_props_.effectiveMobilities((*c_)[cell], (*cmax_)[cell], &cell_viscosity_[np*cell + 0], &cell_relperm_[np*cell + 0], &cell_phasemob_[np*cell + 0]);
        }

        // 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]);
            }
        }
    }

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

//.........这里部分代码省略.........
    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
        phasePressure[phaseIdx].resize(numCells);
    }

    for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
        // we currently hard-code the oil phase as the reference phase!
        assert(pu.phase_used[BlackoilPhases::Liquid]);

        const int opos = pu.phase_pos[BlackoilPhases::Liquid];
        phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];

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

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