C++ ADB::numBlocks方法代码示例

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

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

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

ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
                                 const Cells& cells) const
{
    const int n = cells.size();
    assert(pg.value().size() == n);
    V mu(n);
    V dmudp(n);
    for (int i = 0; i < n; ++i) {
        const double& pg_i = pg.value()[i];
        int regionIdx = cellPvtRegionIdx_[cells[i]];
        double tempInvB = b_[regionIdx](pg_i);
        double tempInvBmu = inverseBmu_[regionIdx](pg_i);
        mu[i] = tempInvB / tempInvBmu;
        dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(pg_i)
                         - tempInvB * inverseBmu_[regionIdx].derivative(pg_i)) / (tempInvBmu * tempInvBmu);
    }

    ADB::M dmudp_diag(dmudp.matrix().asDiagonal());
    const int num_blocks = pg.numBlocks();
    std::vector<ADB::M> jacs(num_blocks);
    for (int block = 0; block < num_blocks; ++block) {
        jacs[block] = dmudp_diag * pg.derivative()[block];
    }
    return ADB::function(std::move(mu), std::move(jacs));
}

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

    /// 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 BlackoilPropsAdFromDeck::bGas(const ADB& pg,
                                      const Cells& cells) const
    {
        if (!phase_usage_.phase_used[Gas]) {
            OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
        }
        const int n = cells.size();
        assert(pg.size() == n);

        V b(n);
        V dbdp(n);
        V dbdr(n);
        const double* rs = 0;

        props_[phase_usage_.phase_pos[Gas]]->b(n, pg.value().data(), rs,
                                               b.data(), dbdp.data(), dbdr.data());

        ADB::M dbdp_diag = spdiag(dbdp);
        const int num_blocks = pg.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dbdp_diag * pg.derivative()[block];
        }
        return ADB::function(b, jacs);
    }

    /// 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]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n formation volume factor values.
    ADB BlackoilPropsAdFromDeck::bOil(const ADB& po,
                                      const ADB& rs,
                                      const Cells& cells) const
    {
        if (!phase_usage_.phase_used[Oil]) {
            OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
        }
        const int n = cells.size();
        assert(po.size() == n);

        V b(n);
        V dbdp(n);
        V dbdr(n);

        props_[phase_usage_.phase_pos[Oil]]->b(n, po.value().data(), rs.value().data(),
                                               b.data(), dbdp.data(), dbdr.data());

        ADB::M dbdp_diag = spdiag(dbdp);
        ADB::M dbdr_diag = spdiag(dbdr);
        const int num_blocks = po.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dbdp_diag * po.derivative()[block] + dbdr_diag * rs.derivative()[block];
        }
        return ADB::function(b, 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
    }

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

 /// 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<ADB> BlackoilPropsAd::relperm(const ADB& sw,
                                           const ADB& so,
                                           const ADB& 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.value().size() == n);
         s_all.col(pu_.phase_pos[Water]) = sw.value();
     }
     if (pu_.phase_used[Oil]) {
         assert(so.value().size() == n);
         s_all.col(pu_.phase_pos[Oil]) = so.value();
     } else {
         OPM_THROW(std::runtime_error, "BlackoilPropsAd::relperm() assumes oil phase is active.");
     }
     if (pu_.phase_used[Gas]) {
         assert(sg.value().size() == n);
         s_all.col(pu_.phase_pos[Gas]) = sg.value();
     }
     Block kr(n, np);
     Block dkr(n, np*np);
     props_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
     const int num_blocks = so.numBlocks();
     std::vector<ADB> relperms;
     relperms.reserve(3);
     typedef const ADB* ADBPtr;
     ADBPtr 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(num_blocks);
             for (int block = 0; block < num_blocks; ++block) {
                 jacs[block] = ADB::M(n, 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 dkr1/ds2.
                 const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
                 ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
                 for (int block = 0; block < num_blocks; ++block) {
                     jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
                 }
             }
             relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
         } else {
             relperms.emplace_back(ADB::null());
         }
     }
     return relperms;
 }

 /// Bubble point curve for Rs as function of oil pressure.
 /// \param[in]  po     Array of n oil pressure values.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n bubble point values for Rs.
 ADB BlackoilPropsAdFromDeck::rsMax(const ADB& po,
                                    const Cells& cells) const
 {
     if (!phase_usage_.phase_used[Oil]) {
         OPM_THROW(std::runtime_error, "Cannot call rsMax(): oil phase not present.");
     }
     const int n = cells.size();
     assert(po.size() == n);
     V rbub(n);
     V drbubdp(n);
     props_[Oil]->rbub(n, po.value().data(), rbub.data(), drbubdp.data());
     ADB::M drbubdp_diag = spdiag(drbubdp);
     const int num_blocks = po.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         jacs[block] = drbubdp_diag * po.derivative()[block];
     }
     return ADB::function(rbub, jacs);
 }

    ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
	                    				      const double* visc) const
    {
	    const int nc = c.size();
    	V inv_mu_w_eff(nc);
    	V dinv_mu_w_eff(nc);
    	for (int i = 0; i < nc; ++i) {
    	    double im = 0, dim = 0;
    	    polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim);
    	    inv_mu_w_eff(i) = im;
    	    dinv_mu_w_eff(i) = dim;
    	}
        ADB::M dim_diag = spdiag(dinv_mu_w_eff);
        const int num_blocks = c.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dim_diag * c.derivative()[block];
        }
        return ADB::function(std::move(inv_mu_w_eff), std::move(jacs));
    }

ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
                                              const Cells& cells,
                                              const std::vector<int>& regionIdx,
                                              const std::vector<NonuniformTableLinear<double>>& tables) const {
    const int n = cells.size();
    assert(X_AD.value().size() == n);
    V x(n);
    V dx(n);
    for (int i = 0; i < n; ++i) {
        const double& X_i = X_AD.value()[i];
        x[i] = tables[regionIdx[cells[i]]](X_i);
        dx[i] = tables[regionIdx[cells[i]]].derivative(X_i);
    }

    ADB::M dx_diag(dx.matrix().asDiagonal());
    const int num_blocks = X_AD.numBlocks();
    std::vector<ADB::M> jacs(num_blocks);
    for (int block = 0; block < num_blocks; ++block) {
        fastSparseProduct(dx_diag, X_AD.derivative()[block], jacs[block]);
    }
    return ADB::function(std::move(x), std::move(jacs));
}

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