Merge pull request #1059 from joakim-hove/create-eclipsegrid

Create eclipsegrid

examples/compute_initial_state.cpp

@@ -92,7 +92,7 @@ try
     Opm::DeckConstPtr deck = parser->parseFile(deck_filename , parseContext);
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
     const double grav = param.getDefault("gravity", unit::gravity);
-    GridManager gm(eclipseState->getInputGrid());
+    GridManager gm(*eclipseState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
     BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
     warnIfUnusedParams(param);

examples/compute_tof.cpp

@@ -177,7 +177,7 @@ try
     EclipseStateConstPtr eclipseState = std::make_shared<EclipseState>(*deck , parseContext);
 
     // Grid init
-    GridManager grid_manager(eclipseState->getInputGrid());
+    GridManager grid_manager(*eclipseState->getInputGrid());
     const UnstructuredGrid& grid = *grid_manager.c_grid();
     // Rock and fluid init
     IncompPropertiesSinglePhase props(deck, eclipseState, grid);

examples/diagnose_relperm.cpp

@@ -78,7 +78,7 @@ try
     Opm::DeckConstPtr deck(parser->parseFile(eclipseFilename, parseContext));
     eclState.reset(new EclipseState(*deck, parseContext));
 
-    GridManager gm(eclState->getInputGrid());
+    GridManager gm(*eclState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
     using boost::filesystem::path; 
     path fpath(eclipseFilename);

examples/wells_example.cpp

@@ -42,7 +42,7 @@ try
     std::cout << "Done!" << std::endl;
 
     // Setup grid
-    GridManager grid(eclipseState->getInputGrid());
+    GridManager grid(*eclipseState->getInputGrid());
 
     // Define rock and fluid properties
     IncompPropertiesFromDeck incomp_properties(deck, eclipseState, *grid.c_grid());

opm/core/grid.h

@@ -242,6 +242,16 @@ struct UnstructuredGrid
        cornerpoint grids.
     */
     int    *cell_facetag;
+
+
+    /*
+       This vector is retained to be able to construct an
+       EclipseGrid representation of the Grid. If the grid
+       processor actually modifies the elements of the zcorn
+       vector from the input the modified version is stored here;
+       otherwise we just use the default.
+    */
+    double * zcorn;
 };
 
 /**
@@ -290,6 +300,15 @@ allocate_grid(size_t ndims     ,
               size_t nnodes    );
 
 
+/**
+   Will allocate storage internally in the grid object to hold a copy
+   of the zcorn data supplied in the second argument.
+*/
+void
+attach_zcorn_copy(struct UnstructuredGrid* G ,
+                  const double * zcorn);
+
+
 /**
  * Import a grid from a character representation stored in file.
  *

opm/core/grid/GridHelpers.cpp

@@ -149,5 +149,29 @@ FaceCellTraits<UnstructuredGrid>::Type faceCells(const UnstructuredGrid& grid)
 {
     return FaceCellsProxy(grid);
 }
+
+
+Opm::EclipseGrid createEclipseGrid(const UnstructuredGrid& grid, const Opm::EclipseGrid& inputGrid ) {
+    const int * dims = UgGridHelpers::cartDims( grid );
+
+    if ((inputGrid.getNX( ) == static_cast<size_t>(dims[0])) &&
+        (inputGrid.getNY( ) == static_cast<size_t>(dims[1])) &&
+        (inputGrid.getNZ( ) == static_cast<size_t>(dims[2]))) {
+        std::vector<int> updatedACTNUM;
+        const int* global_cell = UgGridHelpers::globalCell( grid );
+
+        if (global_cell) {
+            updatedACTNUM.assign( inputGrid.getCartesianSize( ) , 0 );
+            for (int c = 0; c < numCells( grid ); c++) {
+                updatedACTNUM[global_cell[c]] = 1;
+            }
+        }
+
+        return Opm::EclipseGrid( inputGrid, grid.zcorn, updatedACTNUM );
+    } else {
+        throw std::invalid_argument("Size mismatch - dimensions of inputGrid argument and current UnstructuredGrid instance disagree");
+    }
+}
+
 }
 }

opm/core/grid/GridHelpers.hpp

@@ -22,6 +22,7 @@
 #define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
 
 #include <opm/core/grid.h>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
 #include <boost/range/iterator_range.hpp>
@@ -184,6 +185,31 @@ struct CellVolumeIteratorTraits<UnstructuredGrid>
     typedef const double* IteratorType;
 };
 
+    /**
+       Will create an EclipseGrid representation (i.e. based on
+       ZCORN and COORD) of the current UnstructuredGrid
+       instance. When creating the UnstructuredGrid the detailed
+       cornerpoint information is discarded, and it is difficult
+       to go backwards to recreated ZCORN and COORD.
+
+       The current implementation is based on retaining a copy of the
+       zcorn keyword after the Minpvprocessor has modified it.
+
+       We then create a new EclipseGrid instance based on the original
+       input grid, but we "replace" the ZCORN and ACTNUM keywords with
+       the updated versions.
+
+       If the tolerance in the call to create_grid_cornerpoint( ) is
+       finite the grid processing code might collapse cells, the z
+       coordinate transformations from this process will *not* be
+       correctly represented in the EclipseGrid created by this
+       method.
+    */
+
+/// \brief Construct an EclipseGrid instance based on the inputGrid, with modifications to
+/// zcorn and actnum from the dune UnstructuredGrid.
+Opm::EclipseGrid createEclipseGrid(const UnstructuredGrid& grid, const Opm::EclipseGrid& inputGrid );
+
 /// \brief Get an iterator over the cell volumes of a grid positioned at the first cell.
 const double* beginCellVolumes(const UnstructuredGrid& grid);
 

opm/core/grid/GridManager.cpp

@@ -19,6 +19,7 @@
 
 #include "config.h"
 
+#include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/grid/cart_grid.h>
@@ -38,18 +39,18 @@ namespace Opm
 {
 
     /// Construct a 3d corner-point grid from a deck.
-    GridManager::GridManager(Opm::EclipseGridConstPtr eclipseGrid)
+    GridManager::GridManager(const Opm::EclipseGrid& inputGrid)
         : ug_(0)
     {
-        initFromEclipseGrid(eclipseGrid, std::vector<double>());
+        initFromEclipseGrid(inputGrid, std::vector<double>());
     }
 
 
-    GridManager::GridManager(Opm::EclipseGridConstPtr eclipseGrid,
+    GridManager::GridManager(const Opm::EclipseGrid& inputGrid,
                              const std::vector<double>& poreVolumes)
         : ug_(0)
     {
-        initFromEclipseGrid(eclipseGrid, poreVolumes);
+        initFromEclipseGrid(inputGrid, poreVolumes);
     }
 
 
@@ -128,47 +129,49 @@ namespace Opm
 
 
     // Construct corner-point grid from EclipseGrid.
-    void GridManager::initFromEclipseGrid(Opm::EclipseGridConstPtr eclipseGrid,
+    void GridManager::initFromEclipseGrid(const Opm::EclipseGrid& inputGrid,
                                           const std::vector<double>& poreVolumes)
     {
         struct grdecl g;
+        int cells_modified = 0;
         std::vector<int> actnum;
         std::vector<double> coord;
         std::vector<double> zcorn;
         std::vector<double> mapaxes;
 
-        g.dims[0] = eclipseGrid->getNX();
+        g.dims[0] = inputGrid.getNX();
-        g.dims[1] = eclipseGrid->getNY();
+        g.dims[1] = inputGrid.getNY();
-        g.dims[2] = eclipseGrid->getNZ();
+        g.dims[2] = inputGrid.getNZ();
 
-        eclipseGrid->exportMAPAXES( mapaxes );
+        inputGrid.exportMAPAXES( mapaxes );
-        eclipseGrid->exportCOORD( coord );
+        inputGrid.exportCOORD( coord );
-        eclipseGrid->exportZCORN( zcorn );
+        inputGrid.exportZCORN( zcorn );
-        eclipseGrid->exportACTNUM( actnum );
+        inputGrid.exportACTNUM( actnum );
 
         g.coord = coord.data();
         g.zcorn = zcorn.data();
         g.actnum = actnum.data();
         g.mapaxes = mapaxes.data();
 
-        if (!poreVolumes.empty() && (eclipseGrid->getMinpvMode() != MinpvMode::ModeEnum::Inactive)) {
+        if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != MinpvMode::ModeEnum::Inactive)) {
             MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
-            const double minpv_value  = eclipseGrid->getMinpvValue();
+            const double minpv_value  = inputGrid.getMinpvValue();
             // Currently the pinchProcessor is not used and only opmfil is supported
-            //bool opmfil = eclipseGrid->getMinpvMode() == MinpvMode::OpmFIL;
+            //bool opmfil = inputGrid.getMinpvMode() == MinpvMode::OpmFIL;
             bool opmfil = true;
-            mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
+            cells_modified = mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
         }
 
-        const double z_tolerance = eclipseGrid->isPinchActive() ?
+        const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
-            eclipseGrid->getPinchThresholdThickness() : 0.0;
         ug_ = create_grid_cornerpoint(&g, z_tolerance);
         if (!ug_) {
             OPM_THROW(std::runtime_error, "Failed to construct grid.");
         }
+
+        if (cells_modified > 0) {
+            attach_zcorn_copy( ug_ , zcorn.data() );
+        }
     }
-
-
 
 
     void GridManager::createGrdecl(Opm::DeckConstPtr deck, struct grdecl &grdecl)

opm/core/grid/GridManager.hpp

@@ -42,14 +42,14 @@ namespace Opm
     {
     public:
         /// Construct a grid from an EclipseState::EclipseGrid instance.
-        explicit GridManager(Opm::EclipseGridConstPtr eclipseGrid);
+        explicit GridManager(const Opm::EclipseGrid& inputGrid);
 
         /// Construct a grid from an EclipseState::EclipseGrid instance,
         /// giving an explicit set of pore volumes to be used for MINPV
         /// considerations.
         /// \input[in] eclipseGrid    encapsulates a corner-point grid given from a deck
         /// \input[in] poreVolumes    one element per logical cartesian grid element
-        GridManager(Opm::EclipseGridConstPtr eclipseGrid,
+        GridManager(const Opm::EclipseGrid& inputGrid,
                     const std::vector<double>& poreVolumes);
 
         /// Construct a 2d cartesian grid with cells of unit size.
@@ -86,7 +86,7 @@ namespace Opm
         GridManager& operator=(const GridManager& other);
 
         // Construct corner-point grid from EclipseGrid.
-        void initFromEclipseGrid(Opm::EclipseGridConstPtr eclipseGrid,
+        void initFromEclipseGrid(const Opm::EclipseGrid& inputGrid,
                                  const std::vector<double>& poreVolumes);
 
         // The managed UnstructuredGrid.

opm/core/grid/MinpvProcessor.hpp

@@ -48,7 +48,7 @@ namespace Opm
         /// will have the zcorn numbers changed so they are zero-thickness. Any
         /// cell below will be changed to include the deleted volume if mergeMinPCCells is true
         /// els the volume will be lost
-        void process(const std::vector<double>& pv,
+        int process(const std::vector<double>& pv,
                     const double minpv,
                     const std::vector<int>& actnum,
                     const bool mergeMinPVCells,
@@ -61,21 +61,14 @@ namespace Opm
         std::array<int, 3> delta_;
     };
 
-    inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz)
+    inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz) :
-    {
+        dims_( {nx,ny,nz} ),
-        // Not doing init-list init since bracket-init not available
+        delta_( {1 , 2*nx , 4*nx*ny} )
-        // for all compilers we support (gcc 4.4).
+    { }
-        dims_[0] = nx;
-        dims_[1] = ny;
-        dims_[2] = nz;
-        delta_[0] = 1;
-        delta_[1] = 2*nx;
-        delta_[2] = 4*nx*ny;
-    }
 
 
 
-    inline void MinpvProcessor::process(const std::vector<double>& pv,
+    inline int MinpvProcessor::process(const std::vector<double>& pv,
                                         const double minpv,
                                         const std::vector<int>& actnum,
                                         const bool mergeMinPVCells,
@@ -92,6 +85,8 @@ namespace Opm
         //    is true, the higher four zcorn associated with the cell below
         //    is moved to these values (so it gains the deleted volume).
 
+        int cells_modified = 0;
+
         // Check for sane input sizes.
         const size_t log_size = dims_[0] * dims_[1] * dims_[2];
         if (pv.size() != log_size) {
@@ -114,6 +109,7 @@ namespace Opm
                         }
                         setCellZcorn(ii, jj, kk, cz, zcorn);
 
+
                         // optionally add removed volume to the cell below.
                         if (mergeMinPVCells) {
                             // Check if there is a cell below.
@@ -126,10 +122,12 @@ namespace Opm
                                 setCellZcorn(ii, jj, kk + 1, cz_below, zcorn);
                             }
                         }
+                        ++cells_modified;
                     }
                 }
             }
         }
+        return cells_modified;
     }
 
 

opm/core/grid/grid.c

@@ -48,6 +48,8 @@ destroy_grid(struct UnstructuredGrid *g)
 
         free(g->global_cell);
         free(g->cell_facetag);
+
+        free(g->zcorn);
     }
 
     free(g);
@@ -69,6 +71,18 @@ create_grid_empty(void)
 }
 
 
+void
+attach_zcorn_copy(struct UnstructuredGrid* G , const double * zcorn)
+{
+    size_t zcorn_elements = G->cartdims[0] * G->cartdims[1] * G->cartdims[2] * 8;
+    double * new_zcorn = realloc(G->zcorn , zcorn_elements * sizeof * G->zcorn );
+    if (new_zcorn) {
+        G->zcorn = new_zcorn;
+        memcpy(G->zcorn , zcorn , zcorn_elements * sizeof * G->zcorn );
+    }
+}
+
+
 struct UnstructuredGrid *
 allocate_grid(size_t ndims     ,
               size_t ncells    ,
@@ -83,6 +97,9 @@ allocate_grid(size_t ndims     ,
     G = create_grid_empty();
 
     if (G != NULL) {
+        /* zcorn cache - only for output. */
+        G->zcorn = NULL;
+
         /* Grid fields ---------------------------------------- */
         G->dimensions       = ndims;
         G->number_of_cells  = ncells;

tests/test_blackoilstate.cpp

@@ -48,10 +48,10 @@ BOOST_AUTO_TEST_CASE(EqualsDifferentDeckReturnFalse) {
     const auto es2 = Opm::Parser::parse(filename2);
     auto eg2 = es2.getInputGrid();
 
-    GridManager gridManager1(eg1);
+    GridManager gridManager1(*eg1);
     const UnstructuredGrid& grid1 = *gridManager1.c_grid();
 
-    GridManager gridManager2(eg2);
+    GridManager gridManager2(*eg2);
     const UnstructuredGrid& grid2 = *gridManager2.c_grid();
 
     BlackoilState state1( UgGridHelpers::numCells( grid1 ) , UgGridHelpers::numFaces( grid1 ) , 3);
@@ -73,7 +73,7 @@ BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
     std::vector<int> actnum = get_testBlackoilStateActnum();
     eg->resetACTNUM(actnum.data());
 
-    GridManager gridManager(eg);
+    GridManager gridManager(*eg);
     const UnstructuredGrid& grid = *gridManager.c_grid();
 
     BlackoilState state1( UgGridHelpers::numCells( grid ) , UgGridHelpers::numFaces( grid ) , 3);

tests/test_column_extract.cpp

@@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(DisjointColumn)
     for (size_t i = 1; i <= (3 * 3 * 3); i++)
         actnum.push_back(i != 14); // ACTNUM 13*1 0 13* 1
     ep->resetACTNUM(actnum.data());
-    Opm::GridManager manager(ep);
+    Opm::GridManager manager(*ep);
 
     VVI columns;
     Opm::extractColumn(*manager.c_grid(), columns);

tests/test_compressedpropertyaccess.cpp

@@ -60,7 +60,7 @@ struct TestFixture : public Setup
 {
     TestFixture()
         : Setup ()
-        , grid  (ecl->getInputGrid())
+        , grid  (*ecl->getInputGrid())
         , reltol(1.0e-10)
     {
     }

tests/test_pinchprocessor.cpp

@@ -62,7 +62,7 @@ BOOST_AUTO_TEST_CASE(Processing)
 
     const int nc_initial = eclgrid->getNumActive();
 
-    Opm::GridManager gridM(eclgrid, porv);
+    Opm::GridManager gridM(*eclgrid, porv);
     typedef UnstructuredGrid Grid;
     const Grid& grid = *(gridM.c_grid());
     const int* global_cell = Opm::UgGridHelpers::globalCell(grid);

tests/test_relpermdiagnostics.cpp

@@ -58,7 +58,7 @@ BOOST_AUTO_TEST_CASE(diagnosis)
                              });
     Opm::DeckConstPtr deck(parser->parseFile("../tests/relpermDiagnostics.DATA", parseContext));
     eclState.reset(new EclipseState(*deck, parseContext));
-    GridManager gm(eclState->getInputGrid());
+    GridManager gm(*eclState->getInputGrid());
     const UnstructuredGrid& grid = *gm.c_grid();
     std::shared_ptr<CounterLog> counterLog = std::make_shared<CounterLog>(Log::DefaultMessageTypes);
     OpmLog::addBackend( "COUNTERLOG" , counterLog );

tests/test_stoppedwells.cpp

@@ -48,7 +48,7 @@ BOOST_AUTO_TEST_CASE(TestStoppedWells)
     Opm::DeckConstPtr deck(parser->parseFile(filename , parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck , parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     double target_surfacerate_inj;
     double target_surfacerate_prod;

tests/test_ug.cpp

@@ -21,6 +21,7 @@
 #include <opm/core/grid.h>
 #include <opm/core/grid/cornerpoint_grid.h>  /* compute_geometry */
 #include <opm/core/grid/GridManager.hpp>  /* compute_geometry */
+#include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/grid/cpgpreprocess/preprocess.h>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Deck/DeckItem.hpp>
@@ -65,8 +66,8 @@ BOOST_AUTO_TEST_CASE(Equal) {
     BOOST_CHECK( deck1->hasKeyword("ZCORN") );
     BOOST_CHECK( deck1->hasKeyword("COORD") );
     
-    Opm::GridManager grid1(es1.getInputGrid());
+    Opm::GridManager grid1(*es1.getInputGrid());
-    Opm::GridManager grid2(es2.getInputGrid());
+    Opm::GridManager grid2(*es2.getInputGrid());
     
     const UnstructuredGrid* cgrid1 = grid1.c_grid();
     const UnstructuredGrid* cgrid2 = grid2.c_grid();
@@ -89,7 +90,7 @@ BOOST_AUTO_TEST_CASE(EqualEclipseGrid) {
     Opm::EclipseState es(*deck, parseContext);
     auto grid = es.getInputGrid();
 
-    Opm::GridManager gridM(es.getInputGrid());
+    Opm::GridManager gridM(*es.getInputGrid());
     const UnstructuredGrid* cgrid1 = gridM.c_grid();
     struct UnstructuredGrid * cgrid2;
     {
@@ -164,11 +165,13 @@ BOOST_AUTO_TEST_CASE(TOPS_Fully_Specified) {
     Opm::EclipseState es1(*deck1, parseContext);
     Opm::EclipseState es2(*deck2, parseContext);
 
-    Opm::GridManager gridM1(es1.getInputGrid());
+    Opm::GridManager gridM1(*es1.getInputGrid());
-    Opm::GridManager gridM2(es2.getInputGrid());
+    Opm::GridManager gridM2(*es2.getInputGrid());
 
     const UnstructuredGrid* cgrid1 = gridM1.c_grid();
     const UnstructuredGrid* cgrid2 = gridM2.c_grid();
 
     BOOST_CHECK(grid_equal(cgrid1, cgrid2));
+
+    Opm::EclipseGrid grid = Opm::UgGridHelpers::createEclipseGrid( *cgrid1 , *es1.getInputGrid( ) );
 }

tests/test_wellsmanager.cpp

@@ -182,7 +182,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
     Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     {
         Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
@@ -219,7 +219,7 @@ BOOST_AUTO_TEST_CASE(WellsEqual) {
     Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     Opm::WellsManager wellsManager0(eclipseState, 0, *gridManager.c_grid(), NULL);
     Opm::WellsManager wellsManager1(eclipseState, 1, *gridManager.c_grid(), NULL);
@@ -235,7 +235,7 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
     Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     Opm::WellsManager wellsManager0(eclipseState, 0, *gridManager.c_grid(), NULL);
     Opm::WellsManager wellsManager1(eclipseState, 1, *gridManager.c_grid(), NULL);
@@ -258,7 +258,7 @@ BOOST_AUTO_TEST_CASE(WellShutOK) {
     Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     Opm::WellsManager wellsManager2(eclipseState, 2, *gridManager.c_grid(), NULL);
 
@@ -275,7 +275,7 @@ BOOST_AUTO_TEST_CASE(WellSTOPOK) {
     Opm::DeckConstPtr deck(parser->parseFile(filename, parseContext));
 
     Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(*deck, parseContext));
-    Opm::GridManager gridManager(eclipseState->getInputGrid());
+    Opm::GridManager gridManager(*eclipseState->getInputGrid());
 
     Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);