diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c628298f4..2713da0fe 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/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,16 +141,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& 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 pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); - const std::vector 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()); - for (int i = 1; i < max_elems; ++i) { + const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pb.size()); + 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()); - for (int i = 1; i < max_elems; ++i) { + const size_t max_elems = std::min(max_num_cells_faillog, failed_cells_pd.size()); + 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 active_; diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index ae9ab902a..597de98bf 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/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 - void gatherCellDataToGlobalIterator(const Dune::CpGrid& grid, - const ConstIter& local_begin, - const Iter& global_begin) - { -#if HAVE_MPI - FixedSizeIterCopyHandle 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 ))); - } - - 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); - } + const EclipseGrid& inputGrid = eclState().getInputGrid(); + eclIO_.reset(new EclipseIO(eclState(), UgGridHelpers::createEclipseGrid( this->globalGrid() , inputGrid ))); + eclIO_->writeInitial(computeLegacySimProps_(), 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 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( globalSize ), data::TargetType::INIT}; + data::CellData trany = {UnitSystem::measure::transmissibility, std::vector( globalSize ), data::TargetType::INIT}; + data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector( 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(); + const auto& elemEndIt = globalGridView.template end(); + 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(); + const auto& elemEndIt = globalGridView.template end(); + 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_; 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 fluidprops_; - std::unique_ptr geoprops_; std::unique_ptr state_; + NNC nnc_; std::unique_ptr eclIO_; std::unique_ptr output_writer_; boost::any parallel_information_; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 87ff9b3bd..387a055b6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/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> currentFluidInPlace; @@ -419,12 +417,30 @@ protected: const Wells* /* wells */) { } - std::unique_ptr createSolver(const WellModel& well_model) + std::unique_ptr createSolver(WellModel& well_model) { + const auto& gridView = ebosSimulator_.gridView(); + const PhaseUsage& phaseUsage = props_.phaseUsage(); + const std::vector 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(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 legacyCellPvtRegionIdx_; + std::vector legacyPoreVolume_; + std::vector legacyDepth_; typedef RateConverter::SurfaceToReservoirVoidage > 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_;