From 5ebed2f500b668f06328c23e276ded1a8d90500e Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 22 Nov 2016 14:09:06 +0100 Subject: [PATCH 1/6] flow_ebos: only instantiate the grid once it now uses the grid object which gets created by ebos for everything which should make the parallelization efforts easier. I also tried to cut back the use of the legacy property objects (i.e., for the fluid, geologic and rock properties), but this effort ran aground because of the initialization and output code. (also, if those two were fixed, there would probably be issues with the Newton update.) I ran Norne with this and there did not seem to be any notable performance regressions or benefits. --- opm/autodiff/FlowMainEbos.hpp | 157 ++++++++++++++++++++++++++++++++-- opm/autodiff/GridInit.hpp | 30 ++++++- 2 files changed, 177 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index f41fd098a..972193987 100644 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -40,6 +40,7 @@ namespace Opm typedef typename TTAG(EclFlowProblem) TypeTag; typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; // Print startup message if on output rank. void printStartupMessage() @@ -86,8 +87,26 @@ namespace Opm Base::deck_ = ebosSimulator_->gridManager().deck(); Base::eclipse_state_ = ebosSimulator_->gridManager().eclState(); - IOConfig& ioConfig = Base::eclipse_state_->getIOConfig(); - ioConfig.setOutputDir(Base::output_dir_); + + try { + if (Base::output_cout_) { + MissingFeatures::checkKeywords(*Base::deck_); + } + + IOConfig& ioConfig = Base::eclipse_state_->getIOConfig(); + ioConfig.setOutputDir(Base::output_dir_); + + // Possible to force initialization only behavior (NOSIM). + if (Base::param_.has("nosim")) { + const bool nosim = Base::param_.get("nosim"); + ioConfig.overrideNOSIM( nosim ); + } + } + catch (const std::invalid_argument& e) { + std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile_ << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; + throw; + } // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step) if (Base::param_.has("output_interval")) { @@ -95,12 +114,136 @@ namespace Opm eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) ); } + } - // Possible to force initialization only behavior (NOSIM). - if (Base::param_.has("nosim")) { - const bool nosim = Base::param_.get("nosim"); - ioConfig.overrideNOSIM( nosim ); + // Create grid and property objects. + // Writes to: + // grid_init_ + // material_law_manager_ + // fluidprops_ + // rock_comp_ + // gravity_ + // use_local_perm_ + // geoprops_ + void setupGridAndProps() + { + Dune::CpGrid& grid = ebosSimulator_->gridManager().grid(); + Base::material_law_manager_ = ebosSimulator_->problem().materialLawManager(); + grid_init_.reset(new GridInit()); + grid_init_->setGrid(grid); + + // create the legacy properties objects + Base::fluidprops_.reset(new BlackoilPropsAdFromDeck(*deck_, + *eclipse_state_, + Base::material_law_manager_, + grid)); + // Gravity. + assert(UgGridHelpers::dimensions(grid) == 3); + Base::gravity_.fill(0.0); + Base::gravity_[2] = Base::deck_->hasKeyword("NOGRAV") + ? Base::param_.getDefault("gravity", 0.0) + : Base::param_.getDefault("gravity", unit::gravity); + + // Geological properties + Base::use_local_perm_ = true; + geoprops_.reset(new DerivedGeology(grid, + *Base::fluidprops_, + *Base::eclipse_state_, + Base::use_local_perm_, + Base::gravity_.data())); + + } + + // Initialise the reservoir state. Updated fluid props for SWATINIT. + // Writes to: + // state_ + // threshold_pressures_ + // fluidprops_ (if SWATINIT is used) + void setupState() + { + const PhaseUsage pu = Opm::phaseUsageFromDeck(*deck_); + const Grid& grid = Base::grid_init_->grid(); + + // Need old-style fluid object for init purposes (only). + BlackoilPropertiesFromDeck props(*deck_, + *Base::eclipse_state_, + Base::material_law_manager_, + grid.size(/*codim=*/0), + grid.globalCell().data(), + grid.logicalCartesianSize().data(), + param_); + + + // Init state variables (saturation and pressure). + if (param_.has("init_saturation")) { + state_.reset(new ReservoirState(grid.size(/*codim=*/0), + grid.numFaces(), + props.numPhases())); + + initStateBasic(grid.size(/*codim=*/0), + grid.globalCell().data(), + grid.logicalCartesianSize().data(), + grid.numFaces(), + Opm::UgGridHelpers::faceCells(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + Opm::UgGridHelpers::beginCellCentroids(grid), + Grid::dimension, + props, param_, gravity_[2], *state_); + + initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_); + + enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; + if (pu.phase_used[Oil] && pu.phase_used[Gas]) { + const int numPhases = props.numPhases(); + const int numCells = Opm::UgGridHelpers::numCells(grid); + + // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState. + auto& gor = state_->getCellData( BlackoilState::GASOILRATIO ); + const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL ); + for (int c = 0; c < numCells; ++c) { + // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field. + gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]]; + } + } + } else if (deck_->hasKeyword("EQUIL")) { + // Which state class are we really using - what a f... mess? + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); + + initStateEquil(grid, props, *deck_, *Base::eclipse_state_, gravity_[2], *state_); + //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); + } else { + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); + initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::numFaces(grid), + Opm::UgGridHelpers::faceCells(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + Opm::UgGridHelpers::beginCellCentroids(grid), + Opm::UgGridHelpers::dimensions(grid), + props, *deck_, gravity_[2], *state_); } + + // Threshold pressures. + std::map, double> maxDp; + computeMaxDp(maxDp, *deck_, *Base::eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]); + threshold_pressures_ = thresholdPressures(*deck_, *Base::eclipse_state_, grid, maxDp); + std::vector threshold_pressures_nnc = thresholdPressuresNNC(*Base::eclipse_state_, Base::eclipse_state_->getInputNNC(), maxDp); + threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end()); + + // The capillary pressure is scaled in fluidprops_ to match the scaled capillary pressure in props. + if (deck_->hasKeyword("SWATINIT")) { + const int numCells = Opm::UgGridHelpers::numCells(grid); + std::vector cells(numCells); + for (int c = 0; c < numCells; ++c) { cells[c] = c; } + std::vector pc = state_->saturation(); + props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr); + fluidprops_->setSwatInitScaling(state_->saturation(), pc); + } + initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck_->hasKeyword("DISGAS"), deck_->hasKeyword("VAPOIL")); } // Setup linear solver. @@ -123,7 +266,7 @@ namespace Opm Base::param_, *Base::geoprops_, *Base::fluidprops_, - Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr, + /*rockComp=*/nullptr, *Base::fis_solver_, Base::gravity_.data(), Base::deck_->hasKeyword("DISGAS"), diff --git a/opm/autodiff/GridInit.hpp b/opm/autodiff/GridInit.hpp index 5b900ba22..1cbed0a49 100644 --- a/opm/autodiff/GridInit.hpp +++ b/opm/autodiff/GridInit.hpp @@ -73,18 +73,42 @@ namespace Opm class GridInit { public: + GridInit() + { + gridSelfManaged_ = false; + } + /// Initialize from a deck and/or an eclipse state and (logical cartesian) specified pore volumes. GridInit(const EclipseState& eclipse_state, const std::vector& porv) { - grid_.processEclipseFormat(eclipse_state.getInputGrid(), false, false, false, porv); + gridSelfManaged_ = true; + + grid_ = new Dune::CpGrid; + grid_->processEclipseFormat(eclipse_state.getInputGrid(), false, false, false, porv); } + + ~GridInit() + { + if (gridSelfManaged_) + delete grid_; + } + /// Access the created grid. Note that mutable access may be required for load balancing. Dune::CpGrid& grid() { - return grid_; + return *grid_; } + + /// set the grid from the outside + void setGrid(Dune::CpGrid& newGrid) + { + gridSelfManaged_ = false; + grid_ = &newGrid; + } + private: - Dune::CpGrid grid_; + Dune::CpGrid* grid_; + bool gridSelfManaged_; }; #endif // HAVE_OPM_GRID From 1c2a2c417c4406d70529ada0b7b5fe952fbc8648 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Tue, 22 Nov 2016 16:36:23 +0100 Subject: [PATCH 2/6] [bugfix] make flow_ebos work when no wells are present. --- opm/autodiff/BlackoilModelEbos.hpp | 22 +++++++++++++++------- opm/autodiff/StandardWellsDense.hpp | 21 ++++++++++++++------- 2 files changed, 29 insertions(+), 14 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 9bf20a207..f32cf07f4 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -300,7 +300,7 @@ namespace Opm { // Compute the nonlinear update. const int nc = AutoDiffGrid::numCells(grid_); - const int nw = wellModel().wells().number_of_wells; + const int nw = numWells(); BVector x(nc); BVector xw(nw); @@ -447,7 +447,7 @@ namespace Opm { int sizeNonLinear() const { const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int nw = wellModel().wells().number_of_wells; + const int nw = numWells(); return numPhases() * (nc + nw); } @@ -476,8 +476,11 @@ namespace Opm { const auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - // apply well residual to the residual. - wellModel().apply(ebosResid); + if( xw.size() > 0 ) + { + // apply well residual to the residual. + wellModel().apply(ebosResid); + } // set initial guess x = 0.0; @@ -497,9 +500,12 @@ namespace Opm { istlSolver().solve( opA, x, ebosResid ); } - // recover wells. - xw = 0.0; - wellModel().recoverVariable(x, xw); + if( xw.size() > 0 ) + { + // recover wells. + xw = 0.0; + wellModel().recoverVariable(x, xw); + } } //===================================================================== @@ -1249,6 +1255,8 @@ namespace Opm { /// return true if wells are available in the reservoir bool wellsActive() const { return well_model_.wellsActive(); } + int numWells() const { return wellsActive() ? wells().number_of_wells : 0; } + /// return true if wells are available on this process bool localWellsActive() const { return well_model_.localWellsActive(); } diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 0a7951085..c865713b4 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -92,14 +92,17 @@ enum WellVariablePositions { , fluid_(nullptr) , active_(nullptr) , vfp_properties_(nullptr) - , well_perforation_densities_(wells_arg->well_connpos[wells_arg->number_of_wells]) - , well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells]) - , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) - , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) + , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) + , F0_(wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0 ) { - invDuneD_.setBuildMode( Mat::row_wise ); - duneC_.setBuildMode( Mat::row_wise ); - duneB_.setBuildMode( Mat::row_wise ); + if( wells_ ) + { + invDuneD_.setBuildMode( Mat::row_wise ); + duneC_.setBuildMode( Mat::row_wise ); + duneB_.setBuildMode( Mat::row_wise ); + } } void init(const BlackoilPropsAdInterface* fluid_arg, @@ -709,6 +712,10 @@ enum WellVariablePositions { std::vector residual() { + if( ! wellsActive() ) + { + return std::vector(); + } const int np = numPhases(); const int nw = wells().number_of_wells; From c5ca9649d763ab1cedf86a4b3ac8b22a937edbe6 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Tue, 22 Nov 2016 16:36:40 +0100 Subject: [PATCH 3/6] [bugfix] Make initialization work in parallel for flow_ebos. --- opm/autodiff/FlowMain.hpp | 6 +++--- opm/autodiff/FlowMainEbos.hpp | 13 ++++++++++++- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 61bd1f710..2e3d07d7b 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -152,7 +152,7 @@ namespace Opm asImpl().setupOutputWriter(); asImpl().setupLinearSolver(); asImpl().createSimulator(); - + // Run. auto ret = asImpl().runSimulator(); @@ -388,12 +388,12 @@ namespace Opm - // Setup OpmLog backend with output_dir. + // Setup OpmLog backend with output_dir. void setupLogging() { std::string deck_filename = param_.get("deck_filename"); // create logFile - using boost::filesystem::path; + using boost::filesystem::path; path fpath(deck_filename); std::string baseName; std::ostringstream debugFileStream; diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 972193987..b2bfdabb9 100644 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -128,7 +128,18 @@ namespace Opm void setupGridAndProps() { Dune::CpGrid& grid = ebosSimulator_->gridManager().grid(); - Base::material_law_manager_ = ebosSimulator_->problem().materialLawManager(); + + //Base::material_law_manager_ = ebosSimulator_->problem().materialLawManager(); + + grid.switchToGlobalView(); + // Create material law manager. + std::vector compressedToCartesianIdx; + Opm::createGlobalCellArray(grid, compressedToCartesianIdx); + material_law_manager_.reset(new MaterialLawManager()); + material_law_manager_->initFromDeck(*deck_, *eclipse_state_, compressedToCartesianIdx); + + + grid_init_.reset(new GridInit()); grid_init_->setGrid(grid); From e6acf888ccab9114fff06d162ac9172f73927b09 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 30 Nov 2016 13:50:47 +0100 Subject: [PATCH 4/6] flow_ebos: tell the ebos in ourselves to not handle SWATINIT because the flow part also wants to do this. (and it is quite a bit more stubborn!) --- opm/autodiff/BlackoilModelEbos.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index f32cf07f4..3ee5b4db1 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -84,6 +84,10 @@ namespace Properties { NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); SET_BOOL_PROP(EclFlowProblem, DisableWells, true); SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); + +// 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. +SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false); }} namespace Opm { From 2eca5d71e6d8f5b9914f329d835811ca39a6a1f7 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 1 Dec 2016 21:35:37 +0100 Subject: [PATCH 5/6] [bugfix] fix ownerMask for parallel FIP code. --- opm/autodiff/BlackoilModelEbos.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 3ee5b4db1..93d424f29 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -148,7 +148,8 @@ namespace Opm { FIP_PV = 5, //< Pore volume FIP_WEIGHTED_PRESSURE = 6 }; - std::array, 7> fip; + static const int fipValues = FIP_WEIGHTED_PRESSURE + 1 ; + std::array, fipValues> fip; }; // --------- Public methods --------- @@ -1025,7 +1026,7 @@ namespace Opm { const auto& pv = geo_.poreVolume(); const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - for (int i = 0; i<7; i++) { + for (int i = 0; i> values(dims, std::vector(7,0.0)); + std::vector> values(dims, std::vector(FIPData::fipValues,0.0)); std::vector hcpv(dims, 0.0); std::vector pres(dims, 0.0); @@ -1117,11 +1118,12 @@ namespace Opm { // mask[c] is 1 if we need to compute something in parallel const auto & pinfo = boost::any_cast(istlSolver().parallelInformation()); - const auto& mask = pinfo.getOwnerMask(); + const auto& mask = pinfo.updateOwnerMask( fipnum ); + auto comm = pinfo.communicator(); // Compute the global dims value and resize values accordingly. dims = comm.max(dims); - values.resize(dims, std::vector(7,0.0)); + values.resize(dims, std::vector(FIPData::fipValues,0.0)); //Accumulate phases for each region for (int phase = 0; phase < maxnp; ++phase) { From 42ab4d133fd91e49b616771625225bd979f4745a Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 2 Dec 2016 11:13:06 +0100 Subject: [PATCH 6/6] [bugfix] add defunct_well_names to BlackoilModelEbos. --- opm/autodiff/FlowMainEbos.hpp | 3 ++- opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp | 12 ++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index b2bfdabb9..54a6727ce 100644 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -284,7 +284,8 @@ namespace Opm Base::deck_->hasKeyword("VAPOIL"), Base::eclipse_state_, *Base::output_writer_, - Base::threshold_pressures_)); + Base::threshold_pressures_, + Base::defunct_well_names_)); } private: diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index fcb91801c..d4797bc71 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -103,7 +103,8 @@ public: const bool has_vapoil, std::shared_ptr eclipse_state, BlackoilOutputWriterEbos& output_writer, - const std::vector& threshold_pressures_by_face) + const std::vector& threshold_pressures_by_face, + const std::unordered_set& defunct_well_names) : ebosSimulator_(ebosSimulator), param_(param), model_param_(param), @@ -118,6 +119,7 @@ public: terminal_output_(param.getDefault("output_terminal", true)), output_writer_(output_writer), threshold_pressures_by_face_(threshold_pressures_by_face), + defunct_well_names_( defunct_well_names ), is_parallel_run_( false ) { DUNE_UNUSED_PARAMETER(eclipse_state); @@ -229,7 +231,8 @@ public: props_.permeability(), dynamic_list_econ_limited, is_parallel_run_, - well_potentials ); + well_potentials, + defunct_well_names_ ); const Wells* wells = wells_manager.c_wells(); WellState well_state; @@ -733,6 +736,11 @@ protected: std::unique_ptr rateConverter_; // Threshold pressures. std::vector threshold_pressures_by_face_; + // The names of wells that should be defunct + // (e.g. in a parallel run when they are handeled by + // a different process) + std::unordered_set defunct_well_names_; + // Whether this a parallel simulation or not bool is_parallel_run_;