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

    /// Compute the output.
    void CompressibleTpfa::computeResults(BlackoilState& state,
                                          WellState& well_state) const
    {
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        CompletionData completion_data;
        completion_data.wdp = ! wellperf_wdp_.empty() ? const_cast<double*>(&wellperf_wdp_[0]) : 0;
        completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
        completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
        cfs_tpfa_res_wells wells_tmp;
        wells_tmp.W = const_cast<Wells*>(wells_);
        wells_tmp.data = &completion_data;
        cfs_tpfa_res_forces forces;
        forces.wells = &wells_tmp;
        forces.src = NULL;

        double* wpress = ! well_state.bhp      ().empty() ? & well_state.bhp      ()[0] : 0;
        double* wflux  = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;

        cfs_tpfa_res_flux(gg,
                          &forces,
                          props_.numPhases(),
                          &trans_[0],
                          &cell_phasemob_[0],
                          &face_phasemob_[0],
                          &face_gravcap_[0],
                          &state.pressure()[0],
                          wpress,
                          &state.faceflux()[0],
                          wflux);
        cfs_tpfa_res_fpress(gg,
                            props_.numPhases(),
                            &htrans_[0],
                            &face_phasemob_[0],
                            &face_gravcap_[0],
                            h_,
                            &state.pressure()[0],
                            &state.faceflux()[0],
                            &state.facepressure()[0]);

        // Compute well perforation pressures (not done by the C code).
        if (wells_ != 0) {
            const int nw = wells_->number_of_wells;
            for (int w = 0; w < nw; ++w) {
                for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
                    const double bhp = well_state.bhp()[w];
                    well_state.perfPress()[j] = bhp + wellperf_wdp_[j];
                }
            }
        }
    }

    void SimulatorBase<Implementation>::computeWellPotentials(const Wells* wells,
                                                              const BlackoilState& x,
                                                              const WellState& xw,
                                                              std::vector<double>& well_potentials)
    {
        const int nw = wells->number_of_wells;
        const int np = wells->number_of_phases;

        well_potentials.clear();
        well_potentials.resize(nw*np,0.0);       
        for (int w = 0; w < nw; ++w) {
            for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
                const double well_cell_pressure = x.pressure()[wells->well_cells[perf]];
                const double drawdown_used = well_cell_pressure - xw.perfPress()[perf];
                const WellControls* ctrl = wells->ctrls[w];
                const int nwc = well_controls_get_num(ctrl);
                //Loop over all controls until we find a BHP control
                //that specifies what we need...
                double bhp = 0.0;
                for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
                    if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
                        bhp = well_controls_iget_target(ctrl, ctrl_index);
                    }
                    // TODO: do something for thp;
                }
                // Calculate the pressure difference in the well perforation
                const double dp = xw.perfPress()[perf] - xw.bhp()[w];
                const double drawdown_maximum = well_cell_pressure - (bhp + dp);

                for (int phase = 0; phase < np; ++phase) {
                    well_potentials[w*np + phase] += (drawdown_maximum / drawdown_used * xw.perfPhaseRates()[perf*np + phase]);
                }
            }
        }
    }

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

 /// @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 per-solve dynamic properties.
 void CompressibleTpfa::computePerSolveDynamicData(const double /*dt*/,
                                                   const BlackoilState& state,
                                                   const WellState& /*well_state*/)
 {
     computeWellPotentials(state);
     if (rock_comp_props_ && rock_comp_props_->isActive()) {
         computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
     }
 }

            /**
             * Compute average hydrocarbon pressure in all regions.
             *
             * \param[in] state Dynamic reservoir state.
             */
            void
            averagePressure(const BlackoilState& state)
            {
                p_avg_.setZero();

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

                p_avg_ /= ncells_;
            }

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

        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 the residual and Jacobian.
 void CompressibleTpfa::assemble(const double dt,
                                 const BlackoilState& state,
                                 const WellState& well_state)
 {
     const double* cell_press = &state.pressure()[0];
     const double* well_bhp = well_state.bhp().empty() ? NULL : &well_state.bhp()[0];
     const double* z = &state.surfacevol()[0];
     UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
     CompletionData completion_data;
     completion_data.wdp = ! wellperf_wdp_.empty() ? &wellperf_wdp_[0] : 0;
     completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
     completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
     cfs_tpfa_res_wells wells_tmp;
     wells_tmp.W = const_cast<Wells*>(wells_);
     wells_tmp.data = &completion_data;
     cfs_tpfa_res_forces forces;
     forces.wells = &wells_tmp;
     forces.src = NULL; // Check if it is legal to leave it as NULL.
     compr_quantities_gen cq;
     cq.nphases = props_.numPhases();
     cq.Ac = &cell_A_[0];
     cq.dAc = &cell_dA_[0];
     cq.Af = &face_A_[0];
     cq.phasemobf = &face_phasemob_[0];
     cq.voldiscr = &cell_voldisc_[0];
     int was_adjusted = 0;
     if (! (rock_comp_props_ && rock_comp_props_->isActive())) {
         was_adjusted =
             cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                   &face_gravcap_[0], cell_press, well_bhp,
                                   &porevol_[0], h_);
     } else {
         was_adjusted =
             cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                            &face_gravcap_[0], cell_press, well_bhp,
                                            &porevol_[0], &initial_porevol_[0],
                                            &rock_comp_[0], h_);
     }
     singular_ = (was_adjusted == 1);
 }

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

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

    std::cout << "\n================    Test program for weakly compressible two-phase 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");
    EclipseStateConstPtr eclipseState;
    std::unique_ptr<GridManager> grid;
    std::unique_ptr<BlackoilPropertiesInterface> props;
    std::unique_ptr<RockCompressibility> rock_comp;

    ParserPtr parser(new Opm::Parser());
    Opm::DeckConstPtr deck;

    BlackoilState state;
    // bool check_well_controls = false;
    // int max_well_control_iterations = 0;
    double gravity[3] = { 0.0 };
    if (use_deck) {
        ParseMode parseMode;
        std::string deck_filename = param.get<std::string>("deck_filename");
        deck = parser->parseFile(deck_filename , parseMode);
        eclipseState.reset(new EclipseState(deck, parseMode));

        // Grid init
        grid.reset(new GridManager(deck));
        // Rock and fluid init
        props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
        // 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, eclipseState));
        // Gravity.
        gravity[2] = deck->hasKeyword("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);
        } else {
            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
        }
        initBlackoilSurfvol(*grid->c_grid(), *props, state);
    } else {
        // Grid init.
        const int nx = param.getDefault("nx", 100);
        const int ny = param.getDefault("ny", 100);
        const int nz = param.getDefault("nz", 1);
        const double dx = param.getDefault("dx", 1.0);
        const double dy = param.getDefault("dy", 1.0);
        const double dz = param.getDefault("dz", 1.0);
        grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
        // Rock and fluid init.
        props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
        // Rock compressibility.
        rock_comp.reset(new RockCompressibility(param));
        // Gravity.
        gravity[2] = param.getDefault("gravity", 0.0);
        // Init state variables (saturation and pressure).
        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
        initBlackoilSurfvol(*grid->c_grid(), *props, state);
    }

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

    // Initialising src
    int num_cells = grid->c_grid()->number_of_cells;
    std::vector<double> src(num_cells, 0.0);
    if (use_deck) {
        // Do nothing, wells will be the driving force, not source terms.
    } else {
        // Compute pore volumes, in order to enable specifying injection rate
        // terms of total pore volume.
        std::vector<double> porevol;
        if (rock_comp->isActive()) {
            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
        } else {
            computePorevolume(*grid->c_grid(), props->porosity(), porevol);
        }
        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        const double default_injection = use_gravity ? 0.0 : 0.1;
        const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
            *tot_porevol_init/unit::day;
        src[0] = flow_per_sec;
        src[num_cells - 1] = -flow_per_sec;
    }

    // Boundary conditions.
    FlowBCManager bcs;
    if (param.getDefault("use_pside", false)) {
        int pside = param.get<int>("pside");
        double pside_pressure = param.get<double>("pside_pressure");
        bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
    }
//.........这里部分代码省略.........

    /// Compute per-iteration dynamic properties for faces.
    void CompressibleTpfa::computeFaceDynamicData(const double /*dt*/,
                                                  const BlackoilState& state,
                                                  const WellState& /*well_state*/)
    {
        // These are the variables that get computed by this function:
        //
        // std::vector<double> face_A_;
        // std::vector<double> face_phasemob_;
        // std::vector<double> face_gravcap_;
        const int np = props_.numPhases();
        const int nf = grid_.number_of_faces;
        const int dim = grid_.dimensions;
        const double grav = gravity_ ? gravity_[dim - 1] : 0.0;
        std::vector<double> gravcontrib[2];
        std::vector<double> pot[2];
        gravcontrib[0].resize(np);
        gravcontrib[1].resize(np);
        pot[0].resize(np);
        pot[1].resize(np);
        face_A_.resize(nf*np*np);
        face_phasemob_.resize(nf*np);
        face_gravcap_.resize(nf*np);
        for (int face = 0; face < nf; ++face) {
            // Obtain properties from both sides of the face.
            const double face_depth = grid_.face_centroids[face*dim + dim - 1];
            const int* c = &grid_.face_cells[2*face];

            // Get pressures and compute gravity contributions,
            // to decide upwind directions.
            double c_press[2];
            for (int j = 0; j < 2; ++j) {
                if (c[j] >= 0) {
                    // Pressure
                    c_press[j] = state.pressure()[c[j]];
                    // Gravity contribution, gravcontrib = rho*(face_z - cell_z) [per phase].
                    if (grav != 0.0) {
                        const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1];
                        props_.density(1, &cell_A_[np*np*c[j]], &c[j], &gravcontrib[j][0]);
                        for (int p = 0; p < np; ++p) {
                            gravcontrib[j][p] *= depth_diff*grav;
                        }
                    } else {
                        std::fill(gravcontrib[j].begin(), gravcontrib[j].end(), 0.0);
                    }
                } else {
                    // Pressures
                    c_press[j] = state.facepressure()[face];
                    // Gravity contribution.
                    std::fill(gravcontrib[j].begin(), gravcontrib[j].end(), 0.0);
                }
            }

            // Gravity contribution:
            //    gravcapf = rho_1*g*(z_12 - z_1) - rho_2*g*(z_12 - z_2)
            // where _1 and _2 refers to two neigbour cells, z is the
            // z coordinate of the centroid, and z_12 is the face centroid.
            // Also compute the potentials.
            for (int phase = 0; phase < np; ++phase) {
                face_gravcap_[np*face + phase] = gravcontrib[0][phase] - gravcontrib[1][phase];
                pot[0][phase] = c_press[0] + face_gravcap_[np*face + phase];
                pot[1][phase] = c_press[1];
            }

            // Now we can easily find the upwind direction for every phase,
            // we can also tell which boundary faces are inflow bdys.

            // Get upwind mobilities by phase.
            // Get upwind A matrix rows by phase.
            // NOTE:
            // We should be careful to upwind the R factors,
            // the B factors are not that vital.
            //      z = Au = RB^{-1}u,
            // where (this example is for gas-oil)
            //      R = [1 RgL; RoV 1], B = [BL 0 ; 0 BV]
            // (RgL is gas in Liquid phase, RoV is oil in Vapour phase.)
            //      A = [1/BL RgL/BV; RoV/BL 1/BV]
            // This presents us with a dilemma, as V factors should be
            // upwinded according to V phase flow, same for L. What then
            // about the RgL/BV and RoV/BL numbers?
            // We give priority to R, and therefore upwind the rows of A
            // by phase (but remember, Fortran matrix ordering).
            // This prompts the question if we should split the matrix()
            // property method into formation volume and R-factor methods.
            for (int phase = 0; phase < np; ++phase) {
                int upwindc = -1;
                if (c[0] >=0 && c[1] >= 0) {
                    upwindc = (pot[0][phase] < pot[1][phase]) ? c[1] : c[0];
                } else {
                    upwindc = (c[0] >= 0) ? c[0] : c[1];
                }
                face_phasemob_[np*face + phase] = cell_phasemob_[np*upwindc + phase];
                for (int p2 = 0; p2 < np; ++p2) {
                    // Recall: column-major ordering.
                    face_A_[np*np*face + phase + np*p2]
                        = cell_A_[np*np*upwindc + phase + np*p2];
                }
            }
        }
    }

    /// Solve pressure equation, by Newton iterations.
    void CompressibleTpfa::solve(const double dt,
                                 BlackoilState& state,
                                 WellState& well_state)
    {
        const int nc = grid_.number_of_cells;
        const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;

        // Set up dynamic data.
        computePerSolveDynamicData(dt, state, well_state);
        computePerIterationDynamicData(dt, state, well_state);

        // Assemble J and F.
        assemble(dt, state, well_state);

        double inc_norm = 0.0;
        int iter = 0;
        double res_norm = residualNorm();
        std::cout << "\nIteration         Residual        Change in p\n"
                  << std::setw(9) << iter
                  << std::setw(18) << res_norm
                  << std::setw(18) << '*' << std::endl;
        while ((iter < maxiter_) && (res_norm > residual_tol_)) {
            // Solve for increment in Newton method:
            //   incr = x_{n+1} - x_{n} = -J^{-1}F
            // (J is Jacobian matrix, F is residual)
            solveIncrement();
            ++iter;

            // Update pressure vars with increment.
            for (int c = 0; c < nc; ++c) {
                state.pressure()[c] += pressure_increment_[c];
            }
            for (int w = 0; w < nw; ++w) {
                well_state.bhp()[w] += pressure_increment_[nc + w];
            }

            // Stop iterating if increment is small.
            inc_norm = incrementNorm();
            if (inc_norm <= change_tol_) {
                std::cout << std::setw(9) << iter
                          << std::setw(18) << '*'
                          << std::setw(18) << inc_norm << std::endl;
                break;
            }

            // Set up dynamic data.
            computePerIterationDynamicData(dt, state, well_state);

            // Assemble J and F.
            assemble(dt, state, well_state);

            // Update residual norm.
            res_norm = residualNorm();

            std::cout << std::setw(9) << iter
                      << std::setw(18) << res_norm
                      << std::setw(18) << inc_norm << std::endl;
        }

        if ((iter == maxiter_) && (res_norm > residual_tol_) && (inc_norm > change_tol_)) {
            OPM_THROW(std::runtime_error, "CompressibleTpfa::solve() failed to converge in " << maxiter_ << " iterations.");
        }

        std::cout << "Solved pressure in " << iter << " iterations." << std::endl;

        // Compute fluxes and face pressures.
        computeResults(state, well_state);
    }

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