FlowBaseVanguard.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_FLOW_BASE_VANGUARD_HPP
#define OPM_FLOW_BASE_VANGUARD_HPP
 
#include <opm/grid/common/GridEnums.hpp>
#include <opm/grid/common/CartesianIndexMapper.hpp>
#include <opm/grid/LookUpCellCentroid.hh>
 
#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquiferCell.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/io/basevanguard.hh>
#include <opm/models/utils/parametersystem.hh>
#include <opm/models/utils/propertysystem.hh>
 
#include <opm/simulators/flow/BlackoilModelParameters.hpp>
#include <opm/simulators/flow/FlowGenericVanguard.hpp>
 
#include <array>
#include <cstddef>
#include <unordered_map>
#include <vector>
 
namespace Opm {
template <class TypeTag>
class FlowBaseVanguard;
template<typename Grid, typename GridView> struct LookUpCellCentroid;
}
 
namespace Opm::Properties {
 
namespace TTag {
struct FlowBaseVanguard {};
}
 
// declare the properties required by the for the ecl simulator vanguard
template<class TypeTag, class MyTypeTag>
struct EquilGrid {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct EnableOpmRstFile {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct ParsingStrictness {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct SchedRestart {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct EclOutputInterval {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct IgnoreKeywords {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct EdgeWeightsMethod {
    using type = UndefinedProperty;
};
 
#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
template<class TypeTag, class MyTypeTag>
struct NumJacobiBlocks {
    using type = UndefinedProperty;
};
#endif
 
template<class TypeTag, class MyTypeTag>
struct OwnerCellsFirst {
    using type = UndefinedProperty;
};
 
template<class TypeTag, class MyTypeTag>
struct SerialPartitioning {
    using type = UndefinedProperty;
};
 
template<class TypeTag, class MyTypeTag>
struct ZoltanImbalanceTol {
    using type = UndefinedProperty;
};
 
template<class TypeTag, class MyTypeTag>
struct ZoltanParams {
    using type = UndefinedProperty;
};
 
template <class TypeTag, class MyTypeTag>
struct ExternalPartition
{
    using type = UndefinedProperty;
};
 
template<class TypeTag, class MyTypeTag>
struct AllowDistributedWells {
    using type = UndefinedProperty;
};
 
template<class TypeTag>
struct IgnoreKeywords<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr auto value = "";
};
template<class TypeTag>
struct EclDeckFileName<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr auto value = "";
};
template<class TypeTag>
struct EclOutputInterval<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr int value = -1;
};
template<class TypeTag>
struct EnableOpmRstFile<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = false;
};
template<class TypeTag>
struct ParsingStrictness<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr auto value = "normal";
};
template<class TypeTag>
struct SchedRestart<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = false;
};
template<class TypeTag>
struct EdgeWeightsMethod<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr int value = 1;
};
 
#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
template<class TypeTag>
struct NumJacobiBlocks<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr int value = 0;
};
#endif
 
template<class TypeTag>
struct OwnerCellsFirst<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = true;
};
template<class TypeTag>
struct SerialPartitioning<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = false;
};
 
template<class TypeTag>
struct ZoltanImbalanceTol<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr double value = 1.1;
};
 
template<class TypeTag>
struct ZoltanParams<TypeTag,TTag::FlowBaseVanguard> {
    static constexpr auto value = "graph";
};
 
template <class TypeTag>
struct ExternalPartition<TypeTag, TTag::FlowBaseVanguard>
{
    static constexpr auto* value = "";
};
 
template<class TypeTag>
struct AllowDistributedWells<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = false;
};
 
template<class T1, class T2>
struct UseMultisegmentWell;
 
// Same as in BlackoilModelParameters.hpp but for here.
template<class TypeTag>
struct UseMultisegmentWell<TypeTag, TTag::FlowBaseVanguard> {
    static constexpr bool value = true;
};
} // namespace Opm::Properties
 
namespace Opm {
 
template <class TypeTag>
class FlowBaseVanguard : public BaseVanguard<TypeTag>,
                         public FlowGenericVanguard
{
    using ParentType = 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>() };
 
public:
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
 
protected:
    static const int dimension = Grid::dimension;
    static const int dimensionworld = Grid::dimensionworld;
    using Element = typename GridView::template Codim<0>::Entity;
    using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
 
public:
    static void registerParameters()
    {
        Parameters::registerParam<TypeTag, Properties::EclDeckFileName>
            ("The name of the file which contains the ECL deck to be simulated");
        Parameters::registerParam<TypeTag, Properties::EclOutputInterval>
            ("The number of report steps that ought to be skipped between two writes of ECL results");
        Parameters::registerParam<TypeTag, Properties::EnableOpmRstFile>
            ("Include OPM-specific keywords in the ECL restart file to "
             "enable restart of OPM simulators from these files");
        Parameters::registerParam<TypeTag, Properties::IgnoreKeywords>
            ("List of Eclipse keywords which should be ignored. As a ':' separated string.");
        Parameters::registerParam<TypeTag, Properties::ParsingStrictness>
            ("Set strictness of parsing process. Available options are "
             "normal (stop for critical errors), "
             "high (stop for all errors) and "
             "low (as normal, except do not stop due to unsupported "
             "keywords even if marked critical");
        Parameters::registerParam<TypeTag, Properties::SchedRestart>
            ("When restarting: should we try to initialize wells and "
             "groups from historical SCHEDULE section.");
        Parameters::registerParam<TypeTag, Properties::EdgeWeightsMethod>
            ("Choose edge-weighing strategy: 0=uniform, 1=trans, 2=log(trans).");
 
#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
        Parameters::registerParam<TypeTag, Properties::NumJacobiBlocks>
            ("Number of blocks to be created for the Block-Jacobi preconditioner.");
#endif
 
        Parameters::registerParam<TypeTag, Properties::OwnerCellsFirst>
            ("Order cells owned by rank before ghost/overlap cells.");
#if HAVE_MPI
        Parameters::registerParam<TypeTag, Properties::SerialPartitioning>
            ("Perform partitioning for parallel runs on a single process.");
        Parameters::registerParam<TypeTag, Properties::ZoltanImbalanceTol>
            ("Tolerable imbalance of the loadbalancing provided by Zoltan (default: 1.1).");
        Parameters::registerParam<TypeTag, Properties::ZoltanParams>
            ("Configuration of Zoltan partitioner. "
             "Valid options are: graph, hypergraph or scotch. "
             "Alternatively, you can request a configuration to be read "
             "from a JSON file by giving the filename here, ending with '.json.' "
             "See https://sandialabs.github.io/Zoltan/ug_html/ug.html "
             "for available Zoltan options.");
        Parameters::hideParam<TypeTag, Properties::ZoltanParams>();
        Parameters::registerParam<TypeTag, Properties::ExternalPartition>
            ("Name of file from which to load an externally generated "
             "partitioning of the model's active cells for MPI "
             "distribution purposes. If empty, the built-in partitioning "
             "method will be employed.");
        Parameters::hideParam<TypeTag, Properties::ExternalPartition>();
#endif
        Parameters::registerParam<TypeTag, Properties::AllowDistributedWells>
            ("Allow the perforations of a well to be distributed to interior of multiple processes");
        // register here for the use in the tests without BlackoilModelParameters
        Parameters::registerParam<TypeTag, Properties::UseMultisegmentWell>
            ("Use the well model for multi-segment wells instead of the one for single-segment wells");
    }
 
    FlowBaseVanguard(Simulator& simulator)
        : ParentType(simulator)
    {
        fileName_ = Parameters::get<TypeTag, Properties::EclDeckFileName>();
        edgeWeightsMethod_   = Dune::EdgeWeightMethod(Parameters::get<TypeTag, Properties::EdgeWeightsMethod>());
 
#if HAVE_OPENCL || HAVE_ROCSPARSE || HAVE_CUDA
        numJacobiBlocks_ = Parameters::get<TypeTag, Properties::NumJacobiBlocks>();
#endif
 
        ownersFirst_ = Parameters::get<TypeTag, Properties::OwnerCellsFirst>();
#if HAVE_MPI
        serialPartitioning_ = Parameters::get<TypeTag, Properties::SerialPartitioning>();
        zoltanImbalanceTol_ = Parameters::get<TypeTag, Properties::ZoltanImbalanceTol>();
        zoltanParams_ = Parameters::get<TypeTag, Properties::ZoltanParams>();
        externalPartitionFile_ = Parameters::get<TypeTag, Properties::ExternalPartition>();
#endif
        enableDistributedWells_ = Parameters::get<TypeTag, Properties::AllowDistributedWells>();
        ignoredKeywords_ = Parameters::get<TypeTag, Properties::IgnoreKeywords>();
        int output_param = Parameters::get<TypeTag, Properties::EclOutputInterval>();
        if (output_param >= 0)
            outputInterval_ = output_param;
        useMultisegmentWell_ = Parameters::get<TypeTag, Properties::UseMultisegmentWell>();
        enableExperiments_ = enableExperiments;
 
        init();
    }
 
    const CartesianIndexMapper& cartesianMapper() const
    {  return asImp_().cartesianIndexMapper(); }
 
    const std::array<int, dimension>& cartesianDimensions() const
    { return asImp_().cartesianIndexMapper().cartesianDimensions(); }
 
    int cartesianSize() const
    { return asImp_().cartesianIndexMapper().cartesianSize(); }
 
    int equilCartesianSize() const
    { return asImp_().equilCartesianIndexMapper().cartesianSize(); }
 
    unsigned cartesianIndex(unsigned compressedCellIdx) const
    { return asImp_().cartesianIndexMapper().cartesianIndex(compressedCellIdx); }
 
    unsigned cartesianIndex(const std::array<int,dimension>& coords) const
    {
        unsigned cartIndex = coords[0];
        int factor = cartesianDimensions()[0];
        for (unsigned i = 1; i < dimension; ++i) {
            cartIndex += coords[i]*factor;
            factor *= cartesianDimensions()[i];
        }
 
        return cartIndex;
    }
 
    int compressedIndex(int cartesianCellIdx) const
    {
        auto index_pair = cartesianToCompressed_.find(cartesianCellIdx);
        if (index_pair!=cartesianToCompressed_.end())
            return index_pair->second;
        else
            return -1;
    }
 
    int compressedIndexForInterior(int cartesianCellIdx) const
    {
        auto index_pair = cartesianToCompressed_.find(cartesianCellIdx);
        if (index_pair == cartesianToCompressed_.end() ||
            !is_interior_[index_pair->second])
        {
            return -1;
        }
        else
        {
            return index_pair->second;
        }
    }
    void cartesianCoordinate(unsigned cellIdx, std::array<int,3>& ijk) const
    { return asImp_().cartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); }
 
    unsigned equilCartesianIndex(unsigned compressedEquilCellIdx) const
    { return asImp_().equilCartesianIndexMapper().cartesianIndex(compressedEquilCellIdx); }
 
    void equilCartesianCoordinate(unsigned cellIdx, std::array<int,3>& ijk) const
    { return asImp_().equilCartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); }
 
 
    Scalar cellCenterDepth(unsigned globalSpaceIdx) const
    {
        return cellCenterDepth_[globalSpaceIdx];
    }
 
    const std::vector<Scalar>& cellCenterDepths() const
    {
        return cellCenterDepth_;
    }
 
    Scalar cellThickness(unsigned globalSpaceIdx) const
    {
        assert(!cellThickness_.empty());
        return cellThickness_[globalSpaceIdx];
    }
 
    std::size_t globalNumCells() const
    {
        const auto& grid = asImp_().grid();
        if (grid.comm().size() == 1)
        {
            return grid.leafGridView().size(0);
        }
        const auto& gridView = grid.leafGridView();
        constexpr int codim = 0;
        constexpr auto Part = Dune::Interior_Partition;
        auto local_cells = std::distance(gridView.template begin<codim, Part>(),
                                         gridView.template end<codim, Part>());
        return grid.comm().sum(local_cells);
    }
 
    void setupCartesianToCompressed_() {
          this->updateCartesianToCompressedMapping_();
    }
 
protected:
    template<class CartMapper>
    std::function<std::array<double,dimensionworld>(int)>
    cellCentroids_(const CartMapper& cartMapper, const bool& isCpGrid) const
    {
        return [this, cartMapper, isCpGrid](int elemIdx) {
            std::array<double,dimensionworld> centroid;
            const auto rank = this->gridView().comm().rank();
            const auto maxLevel = this->gridView().grid().maxLevel();
            bool useEclipse = !isCpGrid || (isCpGrid && (rank == 0) && (maxLevel == 0));
            if (useEclipse)
            {
                centroid =  this->eclState().getInputGrid().getCellCenter(cartMapper.cartesianIndex(elemIdx));
            }
            else
            {
                LookUpCellCentroid<Grid,GridView> lookUpCellCentroid(this->gridView(), cartMapper, nullptr);
                centroid = lookUpCellCentroid(elemIdx);
            }
            return centroid;
        };
    }
 
    void callImplementationInit()
    {
        asImp_().createGrids_();
        asImp_().filterConnections_();
        std::string outputDir = Parameters::get<TypeTag, Properties::OutputDir>();
        bool enableEclCompatFile = !Parameters::get<TypeTag, Properties::EnableOpmRstFile>();
        asImp_().updateOutputDir_(outputDir, enableEclCompatFile);
        asImp_().finalizeInit_();
    }
 
    void updateCartesianToCompressedMapping_()
    {
        std::size_t num_cells = asImp_().grid().leafGridView().size(0);
        is_interior_.resize(num_cells);
 
        ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
        for (const auto& element : elements(this->gridView()))
        {
            const auto elemIdx = elemMapper.index(element);
            unsigned cartesianCellIdx = cartesianIndex(elemIdx);
            cartesianToCompressed_[cartesianCellIdx] = elemIdx;
            if (element.partitionType() == Dune::InteriorEntity)
            {
                is_interior_[elemIdx] = 1;
            }
            else
            {
                is_interior_[elemIdx] = 0;
            }
        }
    }
 
    void updateCellDepths_()
    {
        int numCells = this->gridView().size(/*codim=*/0);
        cellCenterDepth_.resize(numCells);
 
        ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
 
        const auto num_aqu_cells = this->allAquiferCells();
 
        for(const auto& element : elements(this->gridView())) {
            const unsigned int elemIdx = elemMapper.index(element);
            cellCenterDepth_[elemIdx] = cellCenterDepth(element);
 
            if (!num_aqu_cells.empty()) {
               const unsigned int global_index = cartesianIndex(elemIdx);
               const auto search = num_aqu_cells.find(global_index);
               if (search != num_aqu_cells.end()) {
                    // updating the cell depth using aquifer cell depth
                    cellCenterDepth_[elemIdx] = search->second->depth;
                }
            }
        }
    }
    void updateCellThickness_()
    {
        if (!this->drsdtconEnabled())
            return;
 
        ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
 
        int numElements = this->gridView().size(/*codim=*/0);
        cellThickness_.resize(numElements);
 
        for (const auto& elem : elements(this->gridView())) {
            const unsigned int elemIdx = elemMapper.index(elem);
            cellThickness_[elemIdx] = computeCellThickness(elem);
        }
    }
 
private:
    // computed from averaging cell corner depths
    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);
    }
 
    Scalar computeCellThickness(const typename GridView::template Codim<0>::Entity& element) const
    {
        typedef typename Element::Geometry Geometry;
        static constexpr int zCoord = Element::dimension - 1;
        Scalar zz1 = 0.0;
        Scalar zz2 = 0.0;
 
        const Geometry& geometry = element.geometry();
        // This code only works with CP-grid where the
        // number of corners are 8 and
        // also assumes that the first
        // 4 corners are the top surface and
        // the 4 next are the bottomn.
        assert(geometry.corners() == 8);
        for (int i=0; i < 4; ++i){
            zz1 += geometry.corner(i)[zCoord];
            zz2 += geometry.corner(i+4)[zCoord];
        }
        zz1 /=4;
        zz2 /=4;
        return zz2-zz1;
     }
 
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }
 
    const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }
 
protected:
    std::unordered_map<int,int> cartesianToCompressed_;
 
    std::vector<Scalar> cellCenterDepth_;
 
    std::vector<Scalar> cellThickness_;
 
    std::vector<int> is_interior_;
};
 
} // namespace Opm
 
#endif // OPM_FLOW_BASE_VANGUARD_HPP