mirror of
				https://github.com/OPM/opm-simulators.git
				synced 2025-02-25 18:55:30 -06:00 
			
		
		
		
	Get compressed to cartesian mapping and depths from the vanguard
This commit is contained in:
		@@ -169,6 +169,7 @@ class EclBaseVanguard : public BaseVanguard<TypeTag>
 | 
			
		||||
    using Implementation = GetPropType<TypeTag, Properties::Vanguard>;
 | 
			
		||||
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 | 
			
		||||
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
 | 
			
		||||
    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
 | 
			
		||||
 | 
			
		||||
    enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
 | 
			
		||||
 | 
			
		||||
@@ -178,6 +179,8 @@ public:
 | 
			
		||||
 | 
			
		||||
protected:
 | 
			
		||||
    static const int dimension = Grid::dimension;
 | 
			
		||||
    using Element = typename GridView::template Codim<0>::Entity;
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
public:
 | 
			
		||||
    /*!
 | 
			
		||||
@@ -584,6 +587,16 @@ public:
 | 
			
		||||
        return cartIndex;
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    /*!
 | 
			
		||||
     * \brief Return compressed index from cartesian index
 | 
			
		||||
     *
 | 
			
		||||
     */
 | 
			
		||||
    int compressedIndex(int cartesianCellIdx) const
 | 
			
		||||
    {
 | 
			
		||||
        int index = cartesianToCompressed_[cartesianCellIdx];
 | 
			
		||||
        return index;
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    /*!
 | 
			
		||||
     * \brief Extract Cartesian index triplet (i,j,k) of an active cell.
 | 
			
		||||
     *
 | 
			
		||||
@@ -628,6 +641,17 @@ public:
 | 
			
		||||
        return centroids_;
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    /*!
 | 
			
		||||
     * \brief Returns the depth of an degree of freedom [m]
 | 
			
		||||
     *
 | 
			
		||||
     * For ECL problems this is defined as the average of the depth of an element and is
 | 
			
		||||
     * thus slightly different from the depth of an element's centroid.
 | 
			
		||||
     */
 | 
			
		||||
    Scalar cellCenterDepth(unsigned globalSpaceIdx) const
 | 
			
		||||
    {
 | 
			
		||||
        return cellCenterDepth_[globalSpaceIdx];
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    /*!
 | 
			
		||||
     * \brief Get the number of cells in the global leaf grid view.
 | 
			
		||||
     * \warn This is a collective operation that needs to be called
 | 
			
		||||
@@ -664,6 +688,31 @@ protected:
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
    void updateCartesianToCompressedMapping_()
 | 
			
		||||
    {
 | 
			
		||||
        size_t num_cells = asImp_().grid().leafGridView().size(0);
 | 
			
		||||
        cartesianToCompressed_.resize(cartesianSize(), -1);
 | 
			
		||||
        for (unsigned i = 0; i < num_cells; ++i) {
 | 
			
		||||
            unsigned cartesianCellIdx = cartesianIndex(i);
 | 
			
		||||
            cartesianToCompressed_[cartesianCellIdx] = i;
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    void updateCellDepths_()
 | 
			
		||||
    {
 | 
			
		||||
        int numCells = this->gridView().size(/*codim=*/0);
 | 
			
		||||
        cellCenterDepth_.resize(numCells);
 | 
			
		||||
 | 
			
		||||
        ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
 | 
			
		||||
        auto elemIt = this->gridView().template begin</*codim=*/0>();
 | 
			
		||||
        const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
 | 
			
		||||
        for (; elemIt != elemEndIt; ++elemIt) {
 | 
			
		||||
            const Element& element = *elemIt;
 | 
			
		||||
            const unsigned int elemIdx = elemMapper.index(element);
 | 
			
		||||
            cellCenterDepth_[elemIdx] = cellCenterDepth(element);
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
private:
 | 
			
		||||
    void updateOutputDir_()
 | 
			
		||||
    {
 | 
			
		||||
@@ -693,6 +742,20 @@ private:
 | 
			
		||||
        ioConfig.setEclCompatibleRST(!EWOMS_GET_PARAM(TypeTag, bool, EnableOpmRstFile));
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    Scalar cellCenterDepth(const Element& element) const
 | 
			
		||||
    {
 | 
			
		||||
        typedef typename Element::Geometry Geometry;
 | 
			
		||||
        static constexpr int zCoord = Element::dimension - 1;
 | 
			
		||||
        Scalar zz = 0.0;
 | 
			
		||||
 | 
			
		||||
        const Geometry geometry = element.geometry();
 | 
			
		||||
        const int corners = geometry.corners();
 | 
			
		||||
        for (int i=0; i < corners; ++i)
 | 
			
		||||
            zz += geometry.corner(i)[zCoord];
 | 
			
		||||
 | 
			
		||||
        return zz/Scalar(corners);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    Implementation& asImp_()
 | 
			
		||||
    { return *static_cast<Implementation*>(this); }
 | 
			
		||||
 | 
			
		||||
@@ -734,12 +797,25 @@ protected:
 | 
			
		||||
     * Empty otherwise. Used by EclTransmissibilty.
 | 
			
		||||
     */
 | 
			
		||||
    std::vector<double> centroids_;
 | 
			
		||||
 | 
			
		||||
    /*! \brief Mapping between cartesian and compressed cells.
 | 
			
		||||
     *  It is initialized the first time it is called
 | 
			
		||||
     */
 | 
			
		||||
    std::vector<int> cartesianToCompressed_;
 | 
			
		||||
 | 
			
		||||
    /*! \brief Cell center depths computed
 | 
			
		||||
     *  from averaging cell corner depths
 | 
			
		||||
     */
 | 
			
		||||
    std::vector<Scalar> cellCenterDepth_;
 | 
			
		||||
 | 
			
		||||
    /*! \brief information about wells in parallel
 | 
			
		||||
     *
 | 
			
		||||
     * For each well in the model there is an entry with its name
 | 
			
		||||
     * and a boolean indicating whether it perforates local cells.
 | 
			
		||||
     */
 | 
			
		||||
    std::vector<std::pair<std::string,bool>> parallelWells_;
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <class TypeTag>
 | 
			
		||||
 
 | 
			
		||||
@@ -244,6 +244,10 @@ public:
 | 
			
		||||
 | 
			
		||||
        cartesianIndexMapper_.reset(new CartesianIndexMapper(*grid_));
 | 
			
		||||
        this->updateGridView_();
 | 
			
		||||
        this->updateCartesianToCompressedMapping_();
 | 
			
		||||
        this->updateCellDepths_();
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
#if HAVE_MPI
 | 
			
		||||
        if (mpiSize > 1) {
 | 
			
		||||
            try
 | 
			
		||||
 
 | 
			
		||||
@@ -185,6 +185,8 @@ protected:
 | 
			
		||||
    {
 | 
			
		||||
        grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true));
 | 
			
		||||
        cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
 | 
			
		||||
        this->updateCartesianToCompressedMapping_();
 | 
			
		||||
        this->updateCellDepths_();
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    void filterConnections_()
 | 
			
		||||
 
 | 
			
		||||
@@ -898,7 +898,6 @@ public:
 | 
			
		||||
            maxOilSaturation_.resize(numDof, 0.0);
 | 
			
		||||
        }
 | 
			
		||||
 | 
			
		||||
        updateElementDepths_();
 | 
			
		||||
        readRockParameters_();
 | 
			
		||||
        readMaterialParameters_();
 | 
			
		||||
        readThermalParameters_();
 | 
			
		||||
@@ -1468,9 +1467,10 @@ public:
 | 
			
		||||
    Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
 | 
			
		||||
    {
 | 
			
		||||
        unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
 | 
			
		||||
        return elementCenterDepth_[globalSpaceIdx];
 | 
			
		||||
        return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
    /*!
 | 
			
		||||
     * \copydoc BlackoilProblem::rockCompressibility
 | 
			
		||||
     */
 | 
			
		||||
@@ -2241,40 +2241,6 @@ private:
 | 
			
		||||
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    Scalar cellCenterDepth(const Element& element) const
 | 
			
		||||
    {
 | 
			
		||||
        typedef typename Element::Geometry Geometry;
 | 
			
		||||
        static constexpr int zCoord = Element::dimension - 1;
 | 
			
		||||
        Scalar zz = 0.0;
 | 
			
		||||
 | 
			
		||||
        const Geometry geometry = element.geometry();
 | 
			
		||||
        const int corners = geometry.corners();
 | 
			
		||||
        for (int i=0; i < corners; ++i)
 | 
			
		||||
            zz += geometry.corner(i)[zCoord];
 | 
			
		||||
 | 
			
		||||
        return zz/Scalar(corners);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    void updateElementDepths_()
 | 
			
		||||
    {
 | 
			
		||||
        const auto& simulator = this->simulator();
 | 
			
		||||
        const auto& vanguard = simulator.vanguard();
 | 
			
		||||
        const auto& gridView = vanguard.gridView();
 | 
			
		||||
        const auto& elemMapper = this->elementMapper();;
 | 
			
		||||
 | 
			
		||||
        int numElements = gridView.size(/*codim=*/0);
 | 
			
		||||
        elementCenterDepth_.resize(numElements);
 | 
			
		||||
 | 
			
		||||
        auto elemIt = gridView.template begin</*codim=*/0>();
 | 
			
		||||
        const auto& elemEndIt = gridView.template end</*codim=*/0>();
 | 
			
		||||
        for (; elemIt != elemEndIt; ++elemIt) {
 | 
			
		||||
            const Element& element = *elemIt;
 | 
			
		||||
            const unsigned int elemIdx = elemMapper.index(element);
 | 
			
		||||
 | 
			
		||||
            elementCenterDepth_[elemIdx] = cellCenterDepth(element);
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    // update the parameters needed for DRSDT and DRVDT
 | 
			
		||||
    void updateCompositionChangeLimits_()
 | 
			
		||||
    {
 | 
			
		||||
@@ -2475,7 +2441,7 @@ private:
 | 
			
		||||
                if (!rockTableIdx_.empty()) {
 | 
			
		||||
                    tableIdx = rockTableIdx_[elemIdx];
 | 
			
		||||
                }
 | 
			
		||||
                overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(elementCenterDepth_[elemIdx], /*extrapolation=*/true);
 | 
			
		||||
                overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(vanguard.cellCenterDepth(elemIdx), /*extrapolation=*/true);
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
 | 
			
		||||
@@ -3313,7 +3279,6 @@ private:
 | 
			
		||||
    static std::string briefDescription_;
 | 
			
		||||
 | 
			
		||||
    std::array<std::vector<Scalar>, 2> referencePorosity_;
 | 
			
		||||
    std::vector<Scalar> elementCenterDepth_;
 | 
			
		||||
    EclTransmissibility<TypeTag> transmissibilities_;
 | 
			
		||||
 | 
			
		||||
    std::shared_ptr<EclMaterialLawManager> materialLawManager_;
 | 
			
		||||
 
 | 
			
		||||
@@ -51,10 +51,9 @@ public:
 | 
			
		||||
    using Base::waterCompIdx;
 | 
			
		||||
    using Base::waterPhaseIdx;
 | 
			
		||||
    AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
 | 
			
		||||
                       const std::unordered_map<int, int>& cartesian_to_compressed,
 | 
			
		||||
                       const Simulator& ebosSimulator,
 | 
			
		||||
                       const AquiferCT::AQUCT_data& aquct_data)
 | 
			
		||||
        : Base(aquct_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
 | 
			
		||||
        : Base(aquct_data.aquiferID, connections, ebosSimulator)
 | 
			
		||||
        , aquct_data_(aquct_data)
 | 
			
		||||
    {
 | 
			
		||||
    }
 | 
			
		||||
@@ -92,10 +91,6 @@ protected:
 | 
			
		||||
    // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
 | 
			
		||||
    inline void initializeConnections() override
 | 
			
		||||
    {
 | 
			
		||||
        const auto& eclState = this->ebos_simulator_.vanguard().eclState();
 | 
			
		||||
        const auto& grid = eclState.getInputGrid();
 | 
			
		||||
 | 
			
		||||
        // We hack the cell depth values for now. We can actually get it from elementcontext pos
 | 
			
		||||
        this->cell_depth_.resize(this->size(), this->aquiferDepth());
 | 
			
		||||
        this->alphai_.resize(this->size(), 1.0);
 | 
			
		||||
        this->faceArea_connected_.resize(this->size(), 0.0);
 | 
			
		||||
@@ -108,10 +103,13 @@ protected:
 | 
			
		||||
        this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
 | 
			
		||||
        for (size_t idx = 0; idx < this->size(); ++idx) {
 | 
			
		||||
            const auto global_index = this->connections_[idx].global_index;
 | 
			
		||||
            const int cell_index = this->cartesian_to_compressed_.at(global_index);
 | 
			
		||||
            const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
 | 
			
		||||
 | 
			
		||||
            if (cell_index < 0) //the global_index is not part of this grid
 | 
			
		||||
                continue;
 | 
			
		||||
 | 
			
		||||
            this->cellToConnectionIdx_[cell_index] = idx;
 | 
			
		||||
            const auto cellCenter = grid.getCellCenter(global_index);
 | 
			
		||||
            this->cell_depth_.at(idx) = cellCenter[2];
 | 
			
		||||
            this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
 | 
			
		||||
        }
 | 
			
		||||
        // get default areas for all intersections
 | 
			
		||||
        const auto& gridView = this->ebos_simulator_.vanguard().gridView();
 | 
			
		||||
 
 | 
			
		||||
@@ -53,10 +53,9 @@ public:
 | 
			
		||||
    using Base::waterPhaseIdx;
 | 
			
		||||
 | 
			
		||||
    AquiferFetkovich(const std::vector<Aquancon::AquancCell>& connections,
 | 
			
		||||
                     const std::unordered_map<int, int>& cartesian_to_compressed,
 | 
			
		||||
                     const Simulator& ebosSimulator,
 | 
			
		||||
                     const Aquifetp::AQUFETP_data& aqufetp_data)
 | 
			
		||||
        : Base(aqufetp_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
 | 
			
		||||
        : Base(aqufetp_data.aquiferID, connections, ebosSimulator)
 | 
			
		||||
        , aqufetp_data_(aqufetp_data)
 | 
			
		||||
    {
 | 
			
		||||
    }
 | 
			
		||||
@@ -95,10 +94,6 @@ protected:
 | 
			
		||||
 | 
			
		||||
    inline void initializeConnections() override
 | 
			
		||||
    {
 | 
			
		||||
        const auto& eclState = this->ebos_simulator_.vanguard().eclState();
 | 
			
		||||
        const auto& grid = eclState.getInputGrid();
 | 
			
		||||
 | 
			
		||||
        // We hack the cell depth values for now. We can actually get it from elementcontext pos
 | 
			
		||||
        this->cell_depth_.resize(this->size(), this->aquiferDepth());
 | 
			
		||||
        this->alphai_.resize(this->size(), 1.0);
 | 
			
		||||
        this->faceArea_connected_.resize(this->size(), 0.0);
 | 
			
		||||
@@ -111,10 +106,12 @@ protected:
 | 
			
		||||
        this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
 | 
			
		||||
        for (size_t idx = 0; idx < this->size(); ++idx) {
 | 
			
		||||
            const auto global_index = this->connections_[idx].global_index;
 | 
			
		||||
            const int cell_index = this->cartesian_to_compressed_.at(global_index);
 | 
			
		||||
            const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
 | 
			
		||||
            if (cell_index < 0) //the global_index is not part of this grid
 | 
			
		||||
                continue;
 | 
			
		||||
 | 
			
		||||
            this->cellToConnectionIdx_[cell_index] = idx;
 | 
			
		||||
            const auto cellCenter = grid.getCellCenter(global_index);
 | 
			
		||||
            this->cell_depth_.at(idx) = cellCenter[2];
 | 
			
		||||
            this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
 | 
			
		||||
        }
 | 
			
		||||
        // get areas for all connections
 | 
			
		||||
        const auto& gridView = this->ebos_simulator_.vanguard().gridView();
 | 
			
		||||
 
 | 
			
		||||
@@ -76,12 +76,10 @@ public:
 | 
			
		||||
    // Constructor
 | 
			
		||||
    AquiferInterface(int aqID,
 | 
			
		||||
                     const std::vector<Aquancon::AquancCell>& connections,
 | 
			
		||||
                     const std::unordered_map<int, int>& cartesian_to_compressed,
 | 
			
		||||
                     const Simulator& ebosSimulator)
 | 
			
		||||
        : aquiferID(aqID)
 | 
			
		||||
        , connections_(connections)
 | 
			
		||||
        , ebos_simulator_(ebosSimulator)
 | 
			
		||||
        , cartesian_to_compressed_(cartesian_to_compressed)
 | 
			
		||||
    {
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
@@ -216,7 +214,6 @@ protected:
 | 
			
		||||
    const int aquiferID;
 | 
			
		||||
    const std::vector<Aquancon::AquancCell> connections_;
 | 
			
		||||
    const Simulator& ebos_simulator_;
 | 
			
		||||
    const std::unordered_map<int, int> cartesian_to_compressed_;
 | 
			
		||||
 | 
			
		||||
    // Grid variables
 | 
			
		||||
    std::vector<Scalar> faceArea_connected_;
 | 
			
		||||
 
 | 
			
		||||
@@ -83,7 +83,6 @@ protected:
 | 
			
		||||
 | 
			
		||||
    Simulator& simulator_;
 | 
			
		||||
 | 
			
		||||
    std::unordered_map<int, int> cartesian_to_compressed_;
 | 
			
		||||
    // TODO: probably we can use one variable to store both types of aquifers, because
 | 
			
		||||
    // they share the base class
 | 
			
		||||
    mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy;
 | 
			
		||||
 
 | 
			
		||||
@@ -18,7 +18,6 @@
 | 
			
		||||
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 | 
			
		||||
*/
 | 
			
		||||
 | 
			
		||||
#include <opm/grid/utility/cartesianToCompressed.hpp>
 | 
			
		||||
namespace Opm
 | 
			
		||||
{
 | 
			
		||||
 | 
			
		||||
@@ -170,20 +169,15 @@ BlackoilAquiferModel<TypeTag>::init()
 | 
			
		||||
        throw std::runtime_error("Aquifers currently do not work in parallel.");
 | 
			
		||||
 | 
			
		||||
    // Get all the carter tracy aquifer properties data and put it in aquifers vector
 | 
			
		||||
    const int number_of_cells = simulator_.gridView().size(0);
 | 
			
		||||
 | 
			
		||||
    cartesian_to_compressed_ = cartesianToCompressed(number_of_cells,
 | 
			
		||||
                                                     this->simulator_.vanguard().grid().globalCell().data());
 | 
			
		||||
 | 
			
		||||
    const auto& connections = aquifer.connections();
 | 
			
		||||
    for (const auto& aq : aquifer.ct()) {
 | 
			
		||||
        aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
 | 
			
		||||
                                          cartesian_to_compressed_, this->simulator_, aq);
 | 
			
		||||
                                          this->simulator_, aq);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    for (const auto& aq : aquifer.fetp()) {
 | 
			
		||||
        aquifers_Fetkovich.emplace_back(connections[aq.aquiferID],
 | 
			
		||||
                                        cartesian_to_compressed_, this->simulator_, aq);
 | 
			
		||||
                                        this->simulator_, aq);
 | 
			
		||||
    }
 | 
			
		||||
}
 | 
			
		||||
template <typename TypeTag>
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user