Merge pull request #1019 from andlaus/flow_ebos-remove-geoprops

flow_ebos: remove the legacy geologic properties object

opm/autodiff/BlackoilModelEbos.hpp

@@ -83,6 +83,7 @@ namespace Properties {
 NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
 SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
+SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true);
 
 // SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy
 // code for fluid and satfunc handling gets fully retired.
@@ -90,18 +91,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,7 +141,6 @@ 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
@@ -161,7 +149,6 @@ namespace Opm {
         BlackoilModelEbos(Simulator& ebosSimulator,
                           const ModelParameters& param,
                           const BlackoilPropsAdFromDeck& fluid,
-                          const DerivedGeology&           geo  ,
                           const StandardWellsDense<TypeTag>& well_model,
                           const NewtonIterationBlackoilInterface& linsolver,
                           const bool terminal_output)
@@ -169,7 +156,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 +170,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 +177,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 +863,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 +907,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 +924,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 +1405,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 +1419,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 +1457,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

@@ -56,30 +56,6 @@
 
 namespace Opm
 {
-
-    /// \brief Gather cell data to global random access iterator
-    /// \tparam ConstIter The type of constant iterator.
-    /// \tparam Iter The type of the mutable iterator.
-    /// \param grid The distributed CpGrid where loadbalance has been run.
-    /// \param local The local container from which the data should be sent.
-    /// \param global The global container to gather to.
-    /// \warning The global container has to have the correct size!
-    template<class ConstIter, class Iter>
-    void gatherCellDataToGlobalIterator(const Dune::CpGrid& grid,
-                                         const ConstIter& local_begin,
-                                         const Iter& global_begin)
-    {
-#if HAVE_MPI
-        FixedSizeIterCopyHandle<ConstIter,Iter> handle(local_begin,
-                                                   global_begin);
-        const auto& gatherScatterInf = grid.cellScatterGatherInterface();
-        Dune::VariableSizeCommunicator<> comm(grid.comm(),
-                                              gatherScatterInf);
-        comm.backward(handle);
-#endif
-    }
-
-
     // The FlowMain class is the ebos based black-oil simulator.
     class FlowMainEbos
     {
@@ -87,6 +63,7 @@ namespace Opm
         typedef TTAG(EclFlowProblem) TypeTag;
         typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager;
         typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator;
+        typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
         typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
         typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
         typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
@@ -446,9 +423,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
@@ -598,51 +572,13 @@ namespace Opm
         {
             bool output      = param_.getDefault("output", true);
             bool output_ecl  = param_.getDefault("output_ecl", true);
-            if( output && output_ecl )
+            if( output && output_ecl && grid().comm().rank() == 0 )
             {
-                const Grid& grid = this->globalGrid();
+                exportNncStructure_();
 
-                if( output_cout_ ){
                 const EclipseGrid& inputGrid = eclState().getInputGrid();
-                    eclIO_.reset(new EclipseIO(eclState(), UgGridHelpers::createEclipseGrid( grid , inputGrid )));
+                eclIO_.reset(new EclipseIO(eclState(), UgGridHelpers::createEclipseGrid( this->globalGrid() , inputGrid )));
-                }
+                eclIO_->writeInitial(computeLegacySimProps_(), nnc_);
-
-                const NNC* nnc = &geoprops_->nonCartesianConnections();
-                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());
-                    for( const auto& localkeyval: localtrans)
-                    {
-                        auto& globalval = globaltrans[localkeyval.first].data;
-                        const auto& localval  = localkeyval.second.data;
-
-                        if( output_cout_ )
-                        {
-                            globalval.resize( grid.size(0));
-                        }
-                        gatherCellDataToGlobalIterator(this->grid(), localval.begin(),
-                                                       globalval.begin());
-                    }
-                }
-                else
-                {
-                    globaltrans = geoprops_->simProps(grid);
-                }
-
-                if( output_cout_ )
-                {
-                eclIO_->writeInitial(globaltrans,
-                                              *nnc);
-                }
             }
         }
 
@@ -732,7 +668,6 @@ namespace Opm
             // Create the simulator instance.
             simulator_.reset(new Simulator(*ebosSimulator_,
                                            param_,
-                                           *geoprops_,
                                            *fluidprops_,
                                            *fis_solver_,
                                            FluidSystem::enableDissolvedGas(),
@@ -819,6 +754,153 @@ 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};
+
+            for (size_t i = 0; i < tranx.data.size(); ++i) {
+                tranx.data[0] = 0.0;
+                trany.data[0] = 0.0;
+                tranz.data[0] = 0.0;
+            }
+
+            const Grid& globalGrid = this->globalGrid();
+            const auto& globalGridView = globalGrid.leafGridView();
+            ElementMapper globalElemMapper(globalGridView);
+            const auto& cartesianCellIdx = globalGrid.globalCell();
+
+            const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
+            if (grid().comm().size() < 2) {
+                // in the sequential case we must use the transmissibilites defined by
+                // the problem. (because in the sequential case, the grid manager does
+                // not compute "global" transmissibilities for performance reasons. in
+                // the parallel case, the problem's transmissibilities can't be used
+                // because this object refers to the distributed grid and we need the
+                // sequential version here.)
+                globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
+            }
+
+            auto elemIt = globalGridView.template begin</*codim=*/0>();
+            const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
+            for (; elemIt != elemEndIt; ++ elemIt) {
+                const auto& elem = *elemIt;
+
+                auto isIt = globalGridView.ibegin(elem);
+                const auto& isEndIt = globalGridView.iend(elem);
+                for (; isIt != isEndIt; ++ isIt) {
+                    const auto& is = *isIt;
+
+                    if (!is.neighbor())
+                    {
+                        continue; // intersection is on the domain boundary
+                    }
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+                    unsigned c1 = globalElemMapper.index(is.inside());
+                    unsigned c2 = globalElemMapper.index(is.outside());
+#else
+                    unsigned c1 = globalElemMapper.map(is.inside());
+                    unsigned c2 = globalElemMapper.map(is.outside());
+#endif
+
+                    if (c1 > c2)
+                    {
+                        continue; // we only need to handle each connection once, thank you.
+                    }
+
+
+                    int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]);
+                    int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]);
+                    if (gc2 - gc1 == 1) {
+                        tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
+                    }
+
+                    if (gc2 - gc1 == dims[0]) {
+                        trany.data[gc1] = globalTrans->transmissibility(c1, c2);
+                    }
+
+                    if (gc2 - gc1 == dims[0]*dims[1]) {
+                        tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
+                    }
+                }
+            }
+
+            return {{"TRANX" , tranx},
+                    {"TRANY" , trany} ,
+                    {"TRANZ" , tranz}};
+        }
+
+        void exportNncStructure_()
+        {
+            nnc_ = eclState().getInputNNC();
+            int nx = eclState().getInputGrid().getNX();
+            int ny = eclState().getInputGrid().getNY();
+            //int nz = eclState().getInputGrid().getNZ()
+
+            const Grid& globalGrid = this->globalGrid();
+            const auto& globalGridView = globalGrid.leafGridView();
+            ElementMapper globalElemMapper(globalGridView);
+
+            const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
+            if (grid().comm().size() < 2) {
+                // in the sequential case we must use the transmissibilites defined by
+                // the problem. (because in the sequential case, the grid manager does
+                // not compute "global" transmissibilities for performance reasons. in
+                // the parallel case, the problem's transmissibilities can't be used
+                // because this object refers to the distributed grid and we need the
+                // sequential version here.)
+                globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
+            }
+
+            auto elemIt = globalGridView.template begin</*codim=*/0>();
+            const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
+            for (; elemIt != elemEndIt; ++ elemIt) {
+                const auto& elem = *elemIt;
+
+                auto isIt = globalGridView.ibegin(elem);
+                const auto& isEndIt = globalGridView.iend(elem);
+                for (; isIt != isEndIt; ++ isIt) {
+                    const auto& is = *isIt;
+
+                    if (!is.neighbor())
+                    {
+                        continue; // intersection is on the domain boundary
+                    }
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+                    unsigned c1 = globalElemMapper.index(is.inside());
+                    unsigned c2 = globalElemMapper.index(is.outside());
+#else
+                    unsigned c1 = globalElemMapper.map(is.inside());
+                    unsigned c2 = globalElemMapper.map(is.outside());
+#endif
+
+                    if (c1 > c2)
+                    {
+                        continue; // we only need to handle each connection once, thank you.
+                    }
+
+                    // TODO (?): use the cartesian index mapper to make this code work
+                    // with grids other than Dune::CpGrid. The problem is that we need
+                    // the a mapper for the sequential grid, not for the distributed one.
+                    int cc1 = globalGrid.globalCell()[c1];
+                    int cc2 = globalGrid.globalCell()[c2];
+
+                    if (std::abs(cc1 - cc2) != 1 &&
+                        std::abs(cc1 - cc2) != nx &&
+                        std::abs(cc1 - cc2) != nx*ny)
+                    {
+                        nnc_.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2));
+                    }
+                }
+            }
+        }
+
         std::unique_ptr<EbosSimulator> ebosSimulator_;
         int  mpi_rank_ = 0;
         bool output_cout_ = false;
@@ -827,8 +909,8 @@ 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_;
+        NNC nnc_;
         std::unique_ptr<EclipseIO> eclIO_;
         std::unique_ptr<OutputWriter> output_writer_;
         boost::any parallel_information_;

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_;