C++ Cells::data方法代码示例

 /// Gas formation volume factor.
 /// \param[in]  pg     Array of n gas pressure values.
 /// \param[in]  rv     Array of n vapor oil/gas ratio
 /// \param[in]  cond   Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 ADB BlackoilPropsAd::bGas(const ADB& pg,
                           const ADB& rv,
                           const std::vector<PhasePresence>& /*cond*/,
                           const Cells& cells) const
 {
     if (!pu_.phase_used[Gas]) {
         OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
     }
     const int n = cells.size();
     assert(pg.value().size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     if (pu_.phase_used[Oil]) {
         // Faking a z with the right ratio:
         //   rv = zo/zg
         z.col(pu_.phase_pos[Oil]) = rv.value();
         z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
     }
     Block matrix(n, np*np);
     Block dmatrix(n, np*np);
     props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
     const int phase_ind = pu_.phase_pos[Gas];
     const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
     ADB::M db_diag = spdiag(dmatrix.col(column));
     const int num_blocks = pg.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         jacs[block] = db_diag * pg.derivative()[block];
     }
     return ADB::function(matrix.col(column), jacs);
 }

    /// Water viscosity.
    /// \param[in]  pw     Array of n water pressure values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muWat(const ADB& pw,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muWat(pw.value(), cells), pw.blockPattern());
#else
        if (!pu_.phase_used[Water]) {
            OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
        }
        const int n = cells.size();
        assert(pw.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pw.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Water]));
        const int num_blocks = pw.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pw.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Water]), jacs);
#endif
    }

    /// Gas viscosity.
    /// \param[in]  pg     Array of n gas pressure values.
    /// \param[in]  rv     Array of n vapor oil/gas ratio
    /// \param[in]  cond   Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muGas(const ADB& pg,
                               const ADB& rv,
                               const std::vector<PhasePresence>& cond,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muGas(pg.value(), rv.value(),cond,cells), pg.blockPattern());
#else
        if (!pu_.phase_used[Gas]) {
            OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
        }
        const int n = cells.size();
        assert(pg.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        if (pu_.phase_used[Oil]) {
            // Faking a z with the right ratio:
            //   rv = zo/zg
            z.col(pu_.phase_pos[Oil]) = rv;
            z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
        }
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
        const int num_blocks = pg.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pg.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
#endif
    }

    /// Gas viscosity.
    /// \param[in]  pg     Array of n gas pressure values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muGas(const ADB& pg,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muGas(pg.value(), cells), pg.blockPattern());
#else
        if (!pu_.phase_used[Gas]) {
            THROW("Cannot call muGas(): gas phase not present.");
        }
        const int n = cells.size();
        ASSERT(pg.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
        const int num_blocks = pg.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pg.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
#endif
    }

 /// Relative permeabilities for all phases.
 /// \param[in]  sw     Array of n water saturation values.
 /// \param[in]  so     Array of n oil saturation values.
 /// \param[in]  sg     Array of n gas saturation values.
 /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
 /// \return            An std::vector with 3 elements, each an array of n relperm values,
 ///                    containing krw, kro, krg. Use PhaseIndex for indexing into the result.
 std::vector<V> BlackoilPropsAd::relperm(const V& sw,
                                         const V& so,
                                         const V& sg,
                                         const Cells& cells) const
 {
     const int n = cells.size();
     const int np = props_.numPhases();
     Block s_all(n, np);
     if (pu_.phase_used[Water]) {
         assert(sw.size() == n);
         s_all.col(pu_.phase_pos[Water]) = sw;
     }
     if (pu_.phase_used[Oil]) {
         assert(so.size() == n);
         s_all.col(pu_.phase_pos[Oil]) = so;
     }
     if (pu_.phase_used[Gas]) {
         assert(sg.size() == n);
         s_all.col(pu_.phase_pos[Gas]) = sg;
     }
     Block kr(n, np);
     props_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
     std::vector<V> relperms;
     relperms.reserve(3);
     for (int phase = 0; phase < 3; ++phase) {
         if (pu_.phase_used[phase]) {
             relperms.emplace_back(kr.col(pu_.phase_pos[phase]));
         } else {
             relperms.emplace_back();
         }
     }
     return relperms;
 }

    /// Oil viscosity.
    /// \param[in]  po     Array of n oil pressure values.
    /// \param[in]  rs     Array of n gas solution factor values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muOil(const ADB& po,
                               const ADB& rs,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muOil(po.value(), rs.value(), cells), po.blockPattern());
#else
        if (!pu_.phase_used[Oil]) {
            THROW("Cannot call muOil(): oil phase not present.");
        }
        const int n = cells.size();
        ASSERT(po.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        if (pu_.phase_used[Gas]) {
            // Faking a z with the right ratio:
            //   rs = zg/zo
            z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
            z.col(pu_.phase_pos[Gas]) = rs.value();
        }
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, po.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Oil]));
        const int num_blocks = po.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            // For now, we deliberately ignore the derivative with respect to rs,
            // since the BlackoilPropertiesInterface class does not evaluate it.
            // We would add to the next line: + dmu_drs_diag * rs.derivative()[block]
            jacs[block] = dmu_diag * po.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Oil]), jacs);
#endif
    }

 /// Oil formation volume factor.
 /// \param[in]  po     Array of n oil pressure values.
 /// \param[in]  rs     Array of n gas solution factor values.
 /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 ADB BlackoilPropsAd::bOil(const ADB& po,
                           const ADB& rs,
                           const std::vector<PhasePresence>& /*cond*/,
                           const Cells& cells) const
 {
     if (!pu_.phase_used[Oil]) {
         OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
     }
     const int n = cells.size();
     assert(po.value().size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     if (pu_.phase_used[Gas]) {
         // Faking a z with the right ratio:
         //   rs = zg/zo
         z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
         z.col(pu_.phase_pos[Gas]) = rs.value();
     }
     Block matrix(n, np*np);
     Block dmatrix(n, np*np);
     props_.matrix(n, po.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
     const int phase_ind = pu_.phase_pos[Oil];
     const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
     ADB::M db_diag = spdiag(dmatrix.col(column));
     const int num_blocks = po.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         // For now, we deliberately ignore the derivative with respect to rs,
         // since the BlackoilPropertiesInterface class does not evaluate it.
         // We would add to the next line: + db_drs_diag * rs.derivative()[block]
         jacs[block] = db_diag * po.derivative()[block];
     }
     return ADB::function(matrix.col(column), jacs);
 }

    std::vector<ADB> BlackoilPropsAd::capPress(const ADB& sw,
                                               const ADB& so,
                                               const ADB& sg,
                                               const Cells& cells) const

    {
        const int numCells = cells.size();
        const int numActivePhases = numPhases();
        const int numBlocks = so.numBlocks();

        Block activeSat(numCells, numActivePhases);
        if (pu_.phase_used[Water]) {
            assert(sw.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Water]) = sw.value();
        }
        if (pu_.phase_used[Oil]) {
            assert(so.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Oil]) = so.value();
        } else {
            OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
        }
        if (pu_.phase_used[Gas]) {
            assert(sg.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Gas]) = sg.value();
        }

        Block pc(numCells, numActivePhases);
        Block dpc(numCells, numActivePhases*numActivePhases);
        props_.capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());

        std::vector<ADB> adbCapPressures;
        adbCapPressures.reserve(3);
        const ADB* s[3] = { &sw, &so, &sg };
        for (int phase1 = 0; phase1 < 3; ++phase1) {
            if (pu_.phase_used[phase1]) {
                const int phase1_pos = pu_.phase_pos[phase1];
                std::vector<ADB::M> jacs(numBlocks);
                for (int block = 0; block < numBlocks; ++block) {
                    jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
                }
                for (int phase2 = 0; phase2 < 3; ++phase2) {
                    if (!pu_.phase_used[phase2])
                        continue;
                    const int phase2_pos = pu_.phase_pos[phase2];
                    // Assemble dpc1/ds2.
                    const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
                    ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
                    for (int block = 0; block < numBlocks; ++block) {
                        jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
                    }
                }
                adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
            } else {
                adbCapPressures.emplace_back(ADB::null());
            }
        }
        return adbCapPressures;
    }

 /// Water viscosity.
 /// \param[in]  pw     Array of n water pressure values.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n viscosity values.
 V BlackoilPropsAd::muWat(const V& pw,
                          const Cells& cells) const
 {
     if (!pu_.phase_used[Water]) {
         THROW("Cannot call muWat(): water phase not present.");
     }
     const int n = cells.size();
     ASSERT(pw.size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     Block mu(n, np);
     props_.viscosity(n, pw.data(), z.data(), cells.data(), mu.data(), 0);
     return mu.col(pu_.phase_pos[Water]);
 }

 /// Water formation volume factor.
 /// \param[in]  pw     Array of n water pressure values.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 V BlackoilPropsAd::bWat(const V& pw,
                         const Cells& cells) const
 {
     if (!pu_.phase_used[Water]) {
         OPM_THROW(std::runtime_error, "Cannot call bWat(): water phase not present.");
     }
     const int n = cells.size();
     assert(pw.size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     Block matrix(n, np*np);
     props_.matrix(n, pw.data(), z.data(), cells.data(), matrix.data(), 0);
     const int wi = pu_.phase_pos[Water];
     return matrix.col(wi*np + wi);
 }

 /// Gas formation volume factor.
 /// \param[in]  pg     Array of n gas pressure values.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 ADB BlackoilPropsAd::bGas(const ADB& pg,
                           const Cells& cells) const
 {
     if (!pu_.phase_used[Gas]) {
         THROW("Cannot call muGas(): gas phase not present.");
     }
     const int n = cells.size();
     ASSERT(pg.value().size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     Block matrix(n, np*np);
     Block dmatrix(n, np*np);
     props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
     const int phase_ind = pu_.phase_pos[Gas];
     const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
     ADB::M db_diag = spdiag(dmatrix.col(column));
     const int num_blocks = pg.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         jacs[block] = db_diag * pg.derivative()[block];
     }
     return ADB::function(matrix.col(column), jacs);
 }