Merge pull request #422 from andlaus/satfunc_refactoring

add a few missing "#include <Evaluation.hpp>"

CMakeLists_files.cmake

@@ -114,6 +114,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SolventPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/CPRPreconditioner.hpp
+	opm/autodiff/createGlobalCellArray.hpp
 	opm/autodiff/BlackoilSolventModel.hpp
 	opm/autodiff/BlackoilSolventModel_impl.hpp
 	opm/autodiff/BlackoilSolventState.hpp

examples/flow.cpp

@@ -51,6 +51,7 @@
 #include <opm/core/grid/cornerpoint_grid.h>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/createGlobalCellArray.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/WellsManager.hpp>
@@ -251,15 +252,22 @@ try
     const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
     Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu );
 
+    std::vector<int> compressedToCartesianIdx;
+    Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
+
+    typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
+    auto materialLawManager = std::make_shared<MaterialLawManager>();
+    materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
+
     // Rock and fluid init
-    BlackoilPropertiesFromDeck props( deck, eclipseState,
+    BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager,
                                       Opm::UgGridHelpers::numCells(grid),
                                       Opm::UgGridHelpers::globalCell(grid),
                                       Opm::UgGridHelpers::cartDims(grid),
                                       Opm::UgGridHelpers::beginCellCentroids(grid),
                                       Opm::UgGridHelpers::dimensions(grid), param);
 
-    BlackoilPropsAdFromDeck new_props( deck, eclipseState, grid );
+    BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid );
     // check_well_controls = param.getDefault("check_well_controls", false);
     // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
     // Rock compressibility.

examples/flow_solvent.cpp

@@ -51,6 +51,7 @@
 #include <opm/core/grid/cornerpoint_grid.h>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/createGlobalCellArray.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/WellsManager.hpp>
@@ -247,15 +248,22 @@ try
     const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
     Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu );
 
+    std::vector<int> compressedToCartesianIdx;
+    Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
+
+    typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
+    auto materialLawManager = std::make_shared<MaterialLawManager>();
+    materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
+
     // Rock and fluid init
-    BlackoilPropertiesFromDeck props( deck, eclipseState,
+    BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager,
                                       Opm::UgGridHelpers::numCells(grid),
                                       Opm::UgGridHelpers::globalCell(grid),
                                       Opm::UgGridHelpers::cartDims(grid),
                                       Opm::UgGridHelpers::beginCellCentroids(grid),
                                       Opm::UgGridHelpers::dimensions(grid), param);
 
-    BlackoilPropsAdFromDeck new_props( deck, eclipseState, grid );
+    BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid );
 
     SolventPropsAdFromDeck solvent_props( deck, eclipseState, Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid));
 

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -50,10 +50,11 @@ namespace Opm
     /// Constructor wrapping an opm-core black oil interface.
     BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
                                                      Opm::EclipseStateConstPtr eclState,
+                                                     std::shared_ptr<MaterialLawManager> materialLawManager,
                                                      const UnstructuredGrid& grid,
                                                      const bool init_rock)
     {
-        init(deck, eclState, grid.number_of_cells, grid.global_cell, grid.cartdims, 
+        init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims, 
              grid.cell_centroids, grid.dimensions, init_rock);
     }
 
@@ -64,7 +65,49 @@ namespace Opm
                                                      const Dune::CpGrid& grid,
                                                      const bool init_rock )
     {
-        init(deck, eclState, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
+        auto materialLawManager = std::make_shared<MaterialLawManager>();
+        unsigned number_of_cells = grid.size(0);
+        std::vector<int> compressedToCartesianIdx(number_of_cells);
+        for (unsigned cellIdx = 0; cellIdx < number_of_cells; ++cellIdx) {
+            compressedToCartesianIdx[cellIdx] = grid.globalCell()[cellIdx];
+        }
+        materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
+        init(deck, eclState, materialLawManager, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
+             static_cast<const int*>(&grid.logicalCartesianSize()[0]),
+             grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
+    }
+#endif
+
+    /// Constructor wrapping an opm-core black oil interface.
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
+                                                     Opm::EclipseStateConstPtr eclState,
+                                                     const UnstructuredGrid& grid,
+                                                     const bool init_rock)
+    {
+        auto materialLawManager = std::make_shared<MaterialLawManager>();
+        std::vector<int> compressedToCartesianIdx(grid.number_of_cells);
+        for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) {
+            if (grid.global_cell) {
+                compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx];
+            }
+            else {
+                compressedToCartesianIdx[cellIdx] = cellIdx;
+            }
+        }
+        materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
+        init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims, 
+             grid.cell_centroids, grid.dimensions, init_rock);
+    }
+
+#ifdef HAVE_DUNE_CORNERPOINT
+    /// Constructor wrapping an opm-core black oil interface.
+    BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
+                                                     Opm::EclipseStateConstPtr eclState,
+                                                     std::shared_ptr<MaterialLawManager> materialLawManager,
+                                                     const Dune::CpGrid& grid,
+                                                     const bool init_rock )
+    {
+        init(deck, eclState, materialLawManager, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
              static_cast<const int*>(&grid.logicalCartesianSize()[0]),
              grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
     }
@@ -81,6 +124,9 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
     if (number_of_cells < 0) {
         OPM_THROW(std::runtime_error, "The number of cells is has to be larger than 0.");
     }
+
+    materialLawManager_ = props.materialLawManager_;
+
     // Copy properties that do not depend on the postion within the grid.
     rock_             = props.rock_;
     satprops_         = props.satprops_;
@@ -99,6 +145,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
     template <class CentroidIterator>
     void BlackoilPropsAdFromDeck::init(Opm::DeckConstPtr deck,
                                        Opm::EclipseStateConstPtr eclState,
+                                       std::shared_ptr<MaterialLawManager> materialLawManager,
                                        int number_of_cells,
                                        const int* global_cell,
                                        const int* cart_dims,
@@ -106,6 +153,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
                                        int dimension,
                                        const bool init_rock)
     {
+        materialLawManager_ = materialLawManager;
+
         // retrieve the cell specific PVT table index from the deck
         // and using the grid...
         extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
@@ -239,7 +288,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         SaturationPropsFromDeck* ptr
             = new SaturationPropsFromDeck();
         satprops_.reset(ptr);
-        ptr->init(deck, eclState, number_of_cells, global_cell, begin_cell_centroids, dimension);
+        ptr->init(deck, eclState, materialLawManager_, number_of_cells, global_cell, begin_cell_centroids, dimension);
 
         if (phase_usage_.num_phases != satprops_->numPhases()) {
             OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -59,14 +59,70 @@ namespace Opm
     {
         friend class BlackoilPropsDataHandle;
     public:
-        /// Constructor wrapping an opm-core black oil interface.
+        typedef typename SaturationPropsFromDeck::MaterialLawManager MaterialLawManager;
+
+        /// Constructor to create a blackoil properties from an ECL deck.
+        ///
+        /// The materialLawManager parameter represents the object from opm-material
+        /// which handles the creating and updating parameter objects for the capillary
+        /// pressure/relperm relations for each grid cell. This object is created
+        /// internally for the constructors below, but if it is already available
+        /// externally some performance can be gained by creating it only once.
+        ///
+        /// \param deck The unprocessed ECL deck from opm-parser
+        /// \param eclState The processed ECL deck from opm-parser
+        /// \param materialLawManager The container for the material law parameter objects
+        /// \param grid The grid upon which the simulation is run on.
+        /// \param init_rock If true the rock properties (rock compressibility and
+        ///                  reference pressure) are read from the deck
+        BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
+                                Opm::EclipseStateConstPtr eclState,
+                                std::shared_ptr<MaterialLawManager> materialLawManager,
+                                const UnstructuredGrid& grid,
+                                const bool init_rock = true );
+
+#ifdef HAVE_DUNE_CORNERPOINT
+        /// Constructor to create a blackoil properties from an ECL deck.
+        ///
+        /// The materialLawManager parameter represents the object from opm-material
+        /// which handles the creating and updating parameter objects for the capillary
+        /// pressure/relperm relations for each grid cell. This object is created
+        /// internally for the constructors below, but if it is already available
+        /// externally some performance can be gained by creating it only once.
+        ///
+        /// \param deck The unprocessed ECL deck from opm-parser
+        /// \param eclState The processed ECL deck from opm-parser
+        /// \param materialLawManager The container for the material law parameter objects
+        /// \param grid The grid upon which the simulation is run on.
+        /// \param init_rock If true the rock properties (rock compressibility and
+        ///                  reference pressure) are read from the deck
+        BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
+                                Opm::EclipseStateConstPtr eclState,
+                                std::shared_ptr<MaterialLawManager> materialLawManager,
+                                const Dune::CpGrid& grid,
+                                const bool init_rock = true );
+#endif
+
+        /// Constructor to create a blackoil properties from an ECL deck.
+        ///
+        /// \param deck The unprocessed ECL deck from opm-parser
+        /// \param eclState The processed ECL deck from opm-parser
+        /// \param grid The grid upon which the simulation is run on.
+        /// \param init_rock If true the rock properties (rock compressibility and
+        ///                  reference pressure) are read from the deck
         BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
                                 Opm::EclipseStateConstPtr eclState,
                                 const UnstructuredGrid& grid,
                                 const bool init_rock = true );
 
 #ifdef HAVE_DUNE_CORNERPOINT
-        /// Constructor wrapping an opm-core black oil interface.
+        /// Constructor to create a blackoil properties from an ECL deck.
+        ///
+        /// \param deck The unprocessed ECL deck from opm-parser
+        /// \param eclState The processed ECL deck from opm-parser
+        /// \param grid The grid upon which the simulation is run on.
+        /// \param init_rock If true the rock properties (rock compressibility and
+        ///                  reference pressure) are read from the deck
         BlackoilPropsAdFromDeck(Opm::DeckConstPtr deck,
                                 Opm::EclipseStateConstPtr eclState,
                                 const Dune::CpGrid& grid,
@@ -305,6 +361,7 @@ namespace Opm
         template <class CentroidIterator>
         void init(Opm::DeckConstPtr deck,
                   Opm::EclipseStateConstPtr eclState,
+                  std::shared_ptr<MaterialLawManager> materialLawManager,
                   int number_of_cells,
                   const int* global_cell,
                   const int* cart_dims,
@@ -327,9 +384,11 @@ namespace Opm
         void mapPvtRegions(const std::vector<int>& cells) const;
 
         RockFromDeck rock_;
+
         // This has to be a shared pointer as we must
         // be able to make a copy of *this in the parallel case.
-        std::shared_ptr<SaturationPropsInterface> satprops_;
+        std::shared_ptr<MaterialLawManager> materialLawManager_;
+        std::shared_ptr<SaturationPropsFromDeck> satprops_;
 
         PhaseUsage phase_usage_;
         // bool has_vapoil_;

opm/autodiff/createGlobalCellArray.hpp

@@ -0,0 +1,49 @@
+/*
+  Copyright 2015 Andreas Lauser
+
+  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 3 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/>.
+*/
+#ifndef OPM_CREATE_GLOBAL_CELL_ARRAY_HPP
+#define OPM_CREATE_GLOBAL_CELL_ARRAY_HPP
+
+#include <opm/autodiff/GridHelpers.hpp>
+
+namespace Opm
+{
+/*!
+ * \brief Create a mapping from a global cell index of a grid to the logically
+ *        Cartesian index of the ECL deck.
+ */
+template <class Grid>
+void createGlobalCellArray(const Grid &grid, std::vector<int>& dest)
+{
+    int numCells = Opm::UgGridHelpers::numCells(grid);
+    dest.resize(numCells);
+    const auto& globalCell = Opm::UgGridHelpers::globalCell(grid);
+    std::vector<int> compressedToCartesianIdx(numCells);
+    for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+        if (globalCell) {
+            dest[cellIdx] = globalCell[cellIdx];
+        }
+        else {
+            dest[cellIdx] = cellIdx;
+        }
+    }
+
+}
+}
+
+#endif