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

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

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