flow_ebos: do no longer use the legacy object for geologic properties

it was already almost unused (except for output). Besides making the
overall flow_ebos code leaner because it reduces redundancies, this
patch also implies a small reduduction of memory consumption and a
minor performance improvement. the latter is due to the fact that the
transmissibilities now do not need to be calculated more often than
necessary anymore.

opm/autodiff/BlackoilModelEbos.hpp

@@ -90,18 +90,6 @@ SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false);
 }}
 
 namespace Opm {
-
-
-    class ParameterGroup;
-    class DerivedGeology;
-    class RockCompressibility;
-    class NewtonIterationBlackoilInterface;
-    class VFPProperties;
-    class SimulationDataContainer;
-
-
-
-
     /// A model implementation for three-phase black oil.
     ///
     /// The simulator is capable of handling three-phase problems
@@ -152,16 +140,14 @@ namespace Opm {
         /// \param[in] param            parameters
         /// \param[in] grid             grid data structure
         /// \param[in] fluid            fluid properties
-        /// \param[in] geo              rock properties
         /// \param[in] wells            well structure
         /// \param[in] vfp_properties   Vertical flow performance tables
         /// \param[in] linsolver        linear solver
         /// \param[in] eclState         eclipse state
         /// \param[in] terminal_output  request output to cout/cerr
         BlackoilModelEbos(Simulator& ebosSimulator,
-                          const ModelParameters&          param,
+                          const ModelParameters& param,
                           const BlackoilPropsAdFromDeck& fluid,
-                          const DerivedGeology&           geo  ,
                           const StandardWellsDense<TypeTag>& well_model,
                           const NewtonIterationBlackoilInterface& linsolver,
                           const bool terminal_output)
@@ -169,7 +155,6 @@ namespace Opm {
         , grid_(ebosSimulator_.gridManager().grid())
         , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
         , fluid_ (fluid)
-        , geo_   (geo)
         , vfp_properties_(
             eclState().getTableManager().getVFPInjTables(),
             eclState().getTableManager().getVFPProdTables())
@@ -184,9 +169,6 @@ namespace Opm {
         , dx_old_(AutoDiffGrid::numCells(grid_))
         , isBeginReportStep_(false)
         {
-            const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
-            const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
-            const std::vector<double> depth(geo_.z().data(), geo_.z().data() + geo_.z().size());
             // Wells are active if they are active wells on at least
             // one process.
             int wellsActive = localWellsActive() ? 1 : 0;
@@ -194,7 +176,6 @@ namespace Opm {
             wellModel().setWellsActive( wellsActive );
             // compute global sum of number of cells
             global_nc_ = detail::countGlobalCells(grid_);
-            well_model_.init(fluid_.phaseUsage(), active_, &vfp_properties_, gravity, depth, pv, &rate_converter_, global_nc_);
             if (!istlSolver_)
             {
                 OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
@@ -881,8 +862,6 @@ namespace Opm {
 
             const int nc = Opm::AutoDiffGrid::numCells(grid_);
 
-            const auto& pv = geo_.poreVolume();
-
             const int numComp = numComponents();
             Vector R_sum(numComp);
             Vector B_avg(numComp);
@@ -927,6 +906,16 @@ namespace Opm {
                 }
             }
 
+            Vector pv_vector;
+            pv_vector.resize(nc);
+            const auto& ebosModel = ebosSimulator_.model();
+            const auto& ebosProblem = ebosSimulator_.problem();
+            for (int cellIdx = 0; cellIdx < nc; ++cellIdx) {
+                pv_vector[cellIdx] =
+                    ebosProblem.porosity(cellIdx)*
+                    ebosModel.dofTotalVolume(cellIdx);
+            }
+
             for ( int compIdx = 0; compIdx < numComp; ++compIdx )
             {
                 //tempV.col(compIdx)   = R2.col(compIdx).abs()/pv;
@@ -934,11 +923,10 @@ namespace Opm {
                 Vector& R2_idx    = R2[ compIdx ];
                 for( int cell_idx = 0; cell_idx < nc; ++cell_idx )
                 {
-                    tempV_idx[ cell_idx ] = std::abs( R2_idx[ cell_idx ] ) / pv[ cell_idx ];
+                    tempV_idx[ cell_idx ] = std::abs( R2_idx[ cell_idx ] ) / pv_vector[ cell_idx ];
                 }
             }
 
-            Vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
             Vector wellResidual =  wellModel().residual();
 
             const double pvSum = convergenceReduction(grid_.comm(), global_nc_, numComp,
@@ -1416,8 +1404,8 @@ namespace Opm {
                 std::stringstream errlog;
                 errlog << "Finding the bubble point pressure failed for " << failed_cells_pb.size() << " cells [";
                 errlog << failed_cells_pb[0];
-                const int max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size());
+                const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size());
-                for (int i = 1; i < max_elems; ++i) {
+                for (size_t i = 1; i < max_elems; ++i) {
                     errlog << ", " << failed_cells_pb[i];
                 }
                 if (failed_cells_pb.size() > max_num_cells_faillog) {
@@ -1430,8 +1418,8 @@ namespace Opm {
                 std::stringstream errlog;
                 errlog << "Finding the dew point pressure failed for " << failed_cells_pd.size() << " cells [";
                 errlog << failed_cells_pd[0];
-                const int max_elems = std::min(max_num_cells_faillog, failed_cells_pd.size());
+                const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pd.size());
-                for (int i = 1; i < max_elems; ++i) {
+                for (size_t i = 1; i < max_elems; ++i) {
                     errlog << ", " << failed_cells_pd[i];
                 }
                 if (failed_cells_pd.size() > max_num_cells_faillog) {
@@ -1468,7 +1456,6 @@ namespace Opm {
         const Grid&            grid_;
         const ISTLSolverType*  istlSolver_;
         const BlackoilPropsAdFromDeck& fluid_;
-        const DerivedGeology&           geo_;
         VFPProperties                   vfp_properties_;
         // For each canonical phase -> true if active
         const std::vector<bool>         active_;

opm/autodiff/FlowMainEbos.hpp

@@ -404,9 +404,6 @@ namespace Opm
                                                           materialLawManager(),
                                                           grid));
 
-            // Geological properties
-            bool use_local_perm = param_.getDefault("use_local_perm", true);
-            geoprops_.reset(new DerivedGeology(grid, *fluidprops_, eclState(), use_local_perm, &ebosProblem().gravity()[0]));
         }
 
         const Deck& deck() const
@@ -565,19 +562,16 @@ namespace Opm
                     eclIO_.reset(new EclipseIO(eclState(), UgGridHelpers::createEclipseGrid( grid , inputGrid )));
                 }
 
-                const NNC* nnc = &geoprops_->nonCartesianConnections();
+                eclIO_->writeInitial(computeLegacySimProps_(), eclState().getInputNNC());
+                const NNC* nnc = &eclState().getInputNNC();
                 data::Solution globaltrans;
 
                 if ( must_distribute_ )
                 {
-                    // dirty and dangerous hack!
-                    // We rely on opmfil in GeoProps being hardcoded to true
-                    // which prevents the pinch processing from running.
-                    // Ergo the nncs are unchanged.
                     nnc = &eclState().getInputNNC();
 
                     // Gather the global simProps
-                    data::Solution localtrans = geoprops_->simProps(this->grid());
+                    data::Solution localtrans = computeLegacySimProps_();
                     for( const auto& localkeyval: localtrans)
                     {
                         auto& globalval = globaltrans[localkeyval.first].data;
@@ -593,7 +587,7 @@ namespace Opm
                 }
                 else
                 {
-                    globaltrans = geoprops_->simProps(grid);
+                    globaltrans = computeLegacySimProps_();
                 }
 
                 if( output_cout_ )
@@ -690,7 +684,6 @@ namespace Opm
             // Create the simulator instance.
             simulator_.reset(new Simulator(*ebosSimulator_,
                                            param_,
-                                           *geoprops_,
                                            *fluidprops_,
                                            *fis_solver_,
                                            FluidSystem::enableDissolvedGas(),
@@ -771,6 +764,48 @@ namespace Opm
         std::unordered_set<std::string> defunctWellNames() const
         { return ebosSimulator_->gridManager().defunctWellNames(); }
 
+        data::Solution computeLegacySimProps_()
+        {
+            const int* dims = UgGridHelpers::cartDims(grid());
+            const int globalSize = dims[0]*dims[1]*dims[2];
+
+            data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
+            data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
+            data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
+
+            const auto& eclTrans = ebosSimulator_->problem().eclTransmissibilities();
+            const auto& globalCell = grid().globalCell();
+            size_t num_faces = grid().numFaces();
+            auto fc = UgGridHelpers::faceCells(grid());
+            for (size_t i = 0; i < num_faces; ++i) {
+                auto c1 = std::min(fc(i,0), fc(i,1));
+                auto c2 = std::max(fc(i,0), fc(i,1));
+
+                if (c1 == -1 || c2 == -1)
+                    // boundary
+                    continue;
+
+                int gc1 = globalCell.size()?globalCell[c1]:c1;
+                int gc2 = globalCell.size()?globalCell[c2]:c2;
+
+                if (gc2 - gc1 == 1) {
+                    tranx.data[c1] = eclTrans.transmissibility(c1, c2);
+                }
+
+                if (gc2 - gc1 == dims[0]) {
+                    trany.data[c1] = eclTrans.transmissibility(c1, c2);
+                }
+
+                if (gc2 - gc1 == dims[0]*dims[1]) {
+                    tranz.data[c1] = eclTrans.transmissibility(c1, c2);
+                }
+            }
+            return
+              { {"TRANX" , tranx},
+                {"TRANY" , trany} ,
+                {"TRANZ" , tranz } };
+        }
+
         std::unique_ptr<EbosSimulator> ebosSimulator_;
         int  mpi_rank_ = 0;
         bool output_cout_ = false;
@@ -779,7 +814,6 @@ namespace Opm
         bool output_to_files_ = false;
         std::string output_dir_ = std::string(".");
         std::unique_ptr<BlackoilPropsAdFromDeck> fluidprops_;
-        std::unique_ptr<DerivedGeology> geoprops_;
         std::unique_ptr<ReservoirState> state_;
         std::unique_ptr<EclipseIO> eclIO_;
         std::unique_ptr<OutputWriter> output_writer_;

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -86,7 +86,6 @@ public:
     ///     use_segregation_split (false)  solve for gravity segregation (if false,
     ///                                    segregation is ignored).
     ///
-    /// \param[in] geo           derived geological properties
     /// \param[in] props         fluid and rock properties
     /// \param[in] linsolver     linear solver
     /// \param[in] has_disgas    true for dissolved gas option
@@ -96,7 +95,6 @@ public:
     /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
     SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
                                        const ParameterGroup& param,
-                                       DerivedGeology& geo,
                                        BlackoilPropsAdFromDeck& props,
                                        NewtonIterationBlackoilInterface& linsolver,
                                        const bool has_disgas,
@@ -109,7 +107,6 @@ public:
           model_param_(param),
           solver_param_(param),
           props_(props),
-          geo_(geo),
           solver_(linsolver),
           has_disgas_(has_disgas),
           has_vapoil_(has_vapoil),
@@ -150,6 +147,8 @@ public:
         ExtraData extra;
 
         failureReport_ = SimulatorReport();
+        extractLegacyPoreVolume_();
+        extractLegacyDepth_();
 
         if (output_writer_.isRestart()) {
             // This is a restart, populate WellState and ReservoirState state objects from restart file
@@ -253,8 +252,7 @@ public:
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
 
-            const WellModel well_model(wells, &(wells_manager.wellCollection()), model_param_, terminal_output_);
+            WellModel well_model(wells, &(wells_manager.wellCollection()), model_param_, terminal_output_);
-
             auto solver = createSolver(well_model);
 
             std::vector<std::vector<double>> currentFluidInPlace;
@@ -419,12 +417,30 @@ protected:
                                     const Wells* /* wells */)
     { }
 
-    std::unique_ptr<Solver> createSolver(const WellModel& well_model)
+    std::unique_ptr<Solver> createSolver(WellModel& well_model)
     {
+        const auto& gridView = ebosSimulator_.gridView();
+        const PhaseUsage& phaseUsage = props_.phaseUsage();
+        const std::vector<bool> activePhases = detail::activePhases(phaseUsage);
+        const double gravity = ebosSimulator_.problem().gravity()[2];
+
+        // calculate the number of elements of the compressed sequential grid. this needs
+        // to be done in two steps because the dune communicator expects a reference as
+        // argument for sum()
+        int globalNumCells = gridView.size(/*codim=*/0);
+        globalNumCells = gridView.comm().sum(globalNumCells);
+
+        well_model.init(phaseUsage,
+                        activePhases,
+                        /*vfpProperties=*/nullptr,
+                        gravity,
+                        legacyDepth_,
+                        legacyPoreVolume_,
+                        rateConverter_.get(),
+                        globalNumCells);
         auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
                                                       model_param_,
                                                       props_,
-                                                      geo_,
                                                       well_model,
                                                       solver_,
                                                       terminal_output_));
@@ -807,11 +823,40 @@ protected:
         }
     }
 
+    void extractLegacyPoreVolume_()
+    {
+        const auto& grid = ebosSimulator_.gridManager().grid();
+        const unsigned numCells = grid.size(/*codim=*/0);
+        const auto& ebosProblem = ebosSimulator_.problem();
+        const auto& ebosModel = ebosSimulator_.model();
+
+        legacyPoreVolume_.resize(numCells);
+        for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+            // todo (?): respect rock compressibility
+            legacyPoreVolume_[cellIdx] =
+                ebosModel.dofTotalVolume(cellIdx)
+                *ebosProblem.porosity(cellIdx);
+        }
+    }
+
+    void extractLegacyDepth_()
+    {
+        const auto& grid = ebosSimulator_.gridManager().grid();
+        const unsigned numCells = grid.size(/*codim=*/0);
+
+        legacyDepth_.resize(numCells);
+        for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+            legacyDepth_[cellIdx] =
+                grid.cellCenterDepth(cellIdx);
+        }
+    }
 
     // Data.
     Simulator& ebosSimulator_;
 
     std::vector<int> legacyCellPvtRegionIdx_;
+    std::vector<double> legacyPoreVolume_;
+    std::vector<double> legacyDepth_;
     typedef RateConverter::SurfaceToReservoirVoidage<FluidSystem, std::vector<int> > RateConverterType;
     typedef typename Solver::SolverParameters SolverParameters;
 
@@ -823,8 +868,6 @@ protected:
 
     // Observed objects.
     BlackoilPropsAdFromDeck& props_;
-    // Solvers
-    DerivedGeology& geo_;
     NewtonIterationBlackoilInterface& solver_;
     // Misc. data
     const bool has_disgas_;