diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 90c5c4e2c..508427ea0 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1744,16 +1744,28 @@ typedef Eigen::Array maxResidualAllowed() || CNV[idx] > maxResidualAllowed() || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx)); + const auto msg = std::string("Too large residual for phase ") + materialName(idx); + if (terminal_output_) { + OpmLog::problem(msg); + } + OPM_THROW_NOLOG(Opm::NumericalProblem, msg); } } if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) { - OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation"); + const auto msg = std::string("NaN or too large residual for well control equation"); + if (terminal_output_) { + OpmLog::problem(msg); + } + OPM_THROW_NOLOG(Opm::NumericalProblem, msg); } return converged; @@ -1809,10 +1821,18 @@ typedef Eigen::Array maxResidualAllowed()) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx)); + const auto msg = std::string("Too large residual for phase ") + materialName(idx); + if (terminal_output_) { + OpmLog::problem(msg); + } + OPM_THROW_NOLOG(Opm::NumericalProblem, msg); } } @@ -2198,10 +2218,12 @@ typedef Eigen::ArraywellModel(); + return pressure_solver_.model().wellModel(); } @@ -265,7 +265,7 @@ namespace Opm { /// Return reservoir simulation data (for output functionality) const SimulatorData& getSimulatorData() const { - return transport_model_->getSimulatorData(); + return transport_solver_.model().getSimulatorData(); } diff --git a/opm/autodiff/Compat.hpp b/opm/autodiff/Compat.hpp index d93d71192..9bb208bc1 100644 --- a/opm/autodiff/Compat.hpp +++ b/opm/autodiff/Compat.hpp @@ -21,10 +21,13 @@ #ifndef OPM_SIMULATORS_COMPAT_HPP #define OPM_SIMULATORS_COMPAT_HPP +#include +#include + #include #include #include -#include +#include #include #include #include @@ -167,14 +170,68 @@ inline void solutionToSim( const data::Solution& sol, +inline void wellsToState( const data::Wells& wells, + PhaseUsage phases, + WellStateFullyImplicitBlackoil& state ) { + using rt = data::Rates::opt; -inline void wellsToState( const data::Wells& wells, WellState& state ) { - state.bhp() = wells.bhp; - state.temperature() = wells.temperature; - state.wellRates() = wells.well_rate; - state.perfPress() = wells.perf_pressure; - state.perfRates() = wells.perf_rate; + const auto np = phases.num_phases; + + std::vector< rt > phs( np ); + if( phases.phase_used[BlackoilPhases::Aqua] ) { + phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; + } + + if( phases.phase_used[BlackoilPhases::Liquid] ) { + phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; + } + + if( phases.phase_used[BlackoilPhases::Vapour] ) { + phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; + } + + for( const auto& wm : state.wellMap() ) { + const auto well_index = wm.second[ 0 ]; + const auto& well = wells.at( wm.first ); + + state.bhp()[ well_index ] = well.bhp; + state.temperature()[ well_index ] = well.temperature; + state.currentControls()[ well_index ] = well.control; + + const auto wellrate_index = well_index * np; + for( size_t i = 0; i < phs.size(); ++i ) { + assert( well.rates.has( phs[ i ] ) ); + state.wellRates()[ wellrate_index + i ] = well.rates.get( phs[ i ] ); + } + + const auto perforation_pressure = []( const data::Completion& comp ) { + return comp.pressure; + }; + + const auto perforation_reservoir_rate = []( const data::Completion& comp ) { + return comp.reservoir_rate; + }; + + std::transform( well.completions.begin(), + well.completions.end(), + state.perfPress().begin() + wm.second[ 1 ], + perforation_pressure ); + + std::transform( well.completions.begin(), + well.completions.end(), + state.perfRates().begin() + wm.second[ 1 ], + perforation_reservoir_rate ); + + int local_comp_index = 0; + for (const data::Completion& comp : well.completions) { + const int global_comp_index = wm.second[1] + local_comp_index; + for (int phase_index = 0; phase_index < np; ++phase_index) { + state.perfPhaseRates()[global_comp_index*np + phase_index] = comp.rates.get(phs[phase_index]); + } + ++local_comp_index; + } + } } } diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index efbac9d8f..4110a6df0 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -141,10 +141,10 @@ namespace Opm return EXIT_FAILURE; } asImpl().setupOutput(); - asImpl().setupLogging(); asImpl().readDeckInput(); - asImpl().setupGridAndProps(); + asImpl().setupLogging(); asImpl().extractMessages(); + asImpl().setupGridAndProps(); asImpl().runDiagnostics(); asImpl().setupState(); asImpl().writeInit(); @@ -424,8 +424,16 @@ namespace Opm OpmLog::addBackend( "STREAMLOG", streamLog); std::shared_ptr debugLog = std::make_shared(debugFile, Log::DefaultMessageTypes, false, output_cout_); OpmLog::addBackend( "DEBUGLOG" , debugLog); + const auto& msgLimits = eclipse_state_->getSchedule().getMessageLimits(); + const std::map limits = {{Log::MessageType::Note, msgLimits.getCommentPrintLimit(0)}, + {Log::MessageType::Info, msgLimits.getMessagePrintLimit(0)}, + {Log::MessageType::Warning, msgLimits.getWarningPrintLimit(0)}, + {Log::MessageType::Error, msgLimits.getErrorPrintLimit(0)}, + {Log::MessageType::Problem, msgLimits.getProblemPrintLimit(0)}, + {Log::MessageType::Bug, msgLimits.getBugPrintLimit(0)}}; + prtLog->setMessageLimiter(std::make_shared()); prtLog->setMessageFormatter(std::make_shared(false)); - streamLog->setMessageLimiter(std::make_shared(10)); + streamLog->setMessageLimiter(std::make_shared(10, limits)); streamLog->setMessageFormatter(std::make_shared(true)); // Read parameters. @@ -541,7 +549,7 @@ namespace Opm fluidprops_.reset(new BlackoilPropsAdFromDeck(*deck_, *eclipse_state_, material_law_manager_, grid)); // Rock compressibility. - rock_comp_.reset(new RockCompressibility(*deck_, *eclipse_state_)); + rock_comp_.reset(new RockCompressibility(*deck_, *eclipse_state_, output_cout_)); // Gravity. assert(UgGridHelpers::dimensions(grid) == 3); diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index 0330e9204..83f855833 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -174,6 +174,7 @@ namespace Opm const boost::any& parallelInformation_arg=boost::any()) : iterations_( 0 ), parallelInformation_(parallelInformation_arg), + isIORank_(isIORank(parallelInformation_arg)), parameters_( param ) { } @@ -487,7 +488,11 @@ namespace Opm // Check for failure of linear solver. if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { - OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver."); + const std::string msg("Convergence failure for linear solver."); + if (isIORank_) { + OpmLog::problem(msg); + } + OPM_THROW_NOLOG(LinearSolverProblem, msg); } // Copy solver output to dx. @@ -509,6 +514,7 @@ namespace Opm protected: mutable int iterations_; boost::any parallelInformation_; + bool isIORank_; NewtonIterationBlackoilInterleavedParameters parameters_; }; // end NewtonIterationBlackoilInterleavedImpl @@ -657,7 +663,9 @@ namespace Opm // Check for failure of linear solver. if (!result.converged) { - OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver in computePressureIncrement()."); + const std::string msg("Convergence failure for linear solver in computePressureIncrement()."); + OpmLog::problem(msg); + OPM_THROW_NOLOG(LinearSolverProblem, msg); } // Copy solver output to dx. diff --git a/opm/autodiff/NewtonIterationUtilities.cpp b/opm/autodiff/NewtonIterationUtilities.cpp index 05bf0dc34..7f374ef7e 100644 --- a/opm/autodiff/NewtonIterationUtilities.cpp +++ b/opm/autodiff/NewtonIterationUtilities.cpp @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -288,5 +289,24 @@ namespace Opm + + /// Return true if this is a serial run, or rank zero on an MPI run. + bool isIORank(const boost::any& parallel_info) + { +#if HAVE_MPI + if (parallel_info.type() == typeid(ParallelISTLInformation)) { + const ParallelISTLInformation& info = + boost::any_cast(parallel_info); + return info.communicator().rank() == 0; + } else { + return true; + } +#else + static_cast(parallel_info); // Suppress unused argument warning. + return true; +#endif + } + + } // namespace Opm diff --git a/opm/autodiff/NewtonIterationUtilities.hpp b/opm/autodiff/NewtonIterationUtilities.hpp index 1e383d477..b27106854 100644 --- a/opm/autodiff/NewtonIterationUtilities.hpp +++ b/opm/autodiff/NewtonIterationUtilities.hpp @@ -21,6 +21,7 @@ #define OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED #include +#include #include namespace Opm @@ -58,6 +59,8 @@ namespace Opm Eigen::SparseMatrix& A, AutoDiffBlock::V& b); + /// Return true if this is a serial run, or rank zero on an MPI run. + bool isIORank(const boost::any& parallel_info); } // namespace Opm diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index 77624075b..3c8e010e6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -374,10 +374,12 @@ namespace Opm if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) { std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; } else { + data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...) + combined_sol.insert(simProps.begin(), simProps.end()); // ... insert "extra" data (KR, VISC, ...) eclWriter_->writeTimeStep(timer.reportStepNum(), substep, timer.simulationTimeElapsed(), - simToSolution( state, phaseUsage_ ), + combined_sol, wellState.report(phaseUsage_)); } } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 2becfdd38..4348b9a26 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -277,7 +277,7 @@ namespace Opm /** \brief Whether this process does write to disk */ bool isIORank () const { - parallelOutput_->isIORank(); + return parallelOutput_->isIORank(); } void restore(SimulatorTimerInterface& timer, @@ -420,7 +420,7 @@ namespace Opm Opm::UgGridHelpers::numCells(grid) ); solutionToSim( restarted.first, phaseusage, simulatorstate ); - wellsToState( restarted.second, wellstate ); + wellsToState( restarted.second, phaseusage, wellstate ); } @@ -691,24 +691,26 @@ namespace Opm data::TargetType::SUMMARY ); } if (liquid_active) { + const ADB::V& oipl = sd.fip[Model::SimulatorData::FIP_LIQUID]; + const ADB::V& oipg = vapour_active ? sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL] : ADB::V(); + const ADB::V& oip = vapour_active ? oipl + oipg : oipl; + //Oil in place (liquid phase only) if (hasFRBKeyword(summaryConfig, "OIPL")) { output.insert("OIPL", Opm::UnitSystem::measure::volume, - adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_LIQUID]), + adbVToDoubleVector(oipl), data::TargetType::SUMMARY ); } //Oil in place (gas phase only) if (hasFRBKeyword(summaryConfig, "OIPG")) { output.insert("OIPG", Opm::UnitSystem::measure::volume, - adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]), + adbVToDoubleVector(oipg), data::TargetType::SUMMARY ); } // Oil in place (in liquid and gas phases) if (hasFRBKeyword(summaryConfig, "OIP")) { - ADB::V oip = sd.fip[Model::SimulatorData::FIP_LIQUID] + - sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]; output.insert("OIP", Opm::UnitSystem::measure::volume, adbVToDoubleVector(oip), @@ -716,24 +718,26 @@ namespace Opm } } if (vapour_active) { + const ADB::V& gipg = sd.fip[Model::SimulatorData::FIP_VAPOUR]; + const ADB::V& gipl = liquid_active ? sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS] : ADB::V(); + const ADB::V& gip = liquid_active ? gipg + gipl : gipg; + // Gas in place (gas phase only) if (hasFRBKeyword(summaryConfig, "GIPG")) { output.insert("GIPG", Opm::UnitSystem::measure::volume, - adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPOUR]), + adbVToDoubleVector(gipg), data::TargetType::SUMMARY ); } // Gas in place (liquid phase only) if (hasFRBKeyword(summaryConfig, "GIPL")) { output.insert("GIPL", Opm::UnitSystem::measure::volume, - adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]), + adbVToDoubleVector(gipl), data::TargetType::SUMMARY ); } // Gas in place (in both liquid and gas phases) if (hasFRBKeyword(summaryConfig, "GIP")) { - ADB::V gip = sd.fip[Model::SimulatorData::FIP_VAPOUR] + - sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]; output.insert("GIP", Opm::UnitSystem::measure::volume, adbVToDoubleVector(gip), diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 67a5951c4..ab22e9061 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -280,9 +280,23 @@ namespace Opm data::Wells res = WellState::report(pu); const int nw = this->numWells(); - // If there are now wells numPhases throws a floating point - // exception. - const int np = nw ? this->numPhases() : -1; + if( nw == 0 ) return res; + const int np = pu.num_phases; + + + using rt = data::Rates::opt; + std::vector< rt > phs( np ); + if( pu.phase_used[BlackoilPhases::Aqua] ) { + phs.at( pu.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; + } + + if( pu.phase_used[BlackoilPhases::Liquid] ) { + phs.at( pu.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; + } + + if( pu.phase_used[BlackoilPhases::Vapour] ) { + phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; + } /* this is a reference or example on **how** to convert from * WellState to something understood by opm-output. it is intended @@ -291,42 +305,23 @@ namespace Opm * representations. */ - for( auto w = 0; w < nw; ++w ) { - using rt = data::Rates::opt; - std::map< size_t, data::Completion > completions; + for( const auto& wt : this->wellMap() ) { + const auto w = wt.second[ 0 ]; + auto& well = res.at( wt.first ); + well.control = this->currentControls()[ w ]; - // completions aren't supported yet - //const auto* begin = wells_->well_connpos + w; - //const auto* end = wells_->well_connpos + w + 1; - //for( auto* i = begin; i != end; ++i ) { - // const auto perfrate = this->perfPhaseRates().begin() + *i; - // data::Rates perfrates; - // perfrates.set( rt::wat, *(perfrate + 0) ); - // perfrates.set( rt::oil, *(perfrate + 1) ); - // perfrates.set( rt::gas, *(perfrate + 2) ); + int local_comp_index = 0; + for( auto& comp : well.completions ) { + const auto rates = this->perfPhaseRates().begin() + + (np * wt.second[ 1 ]) + + (np * local_comp_index); + ++local_comp_index; - // const size_t active_index = wells_->well_cells[ *i ]; - - // completions.emplace( active_index, - // data::Completion{ active_index, perfrates } ); - //} - - const auto wellrate_index = np * w; - const auto& wv = this->wellRates(); - - data::Rates wellrates; - if( np == 3 ) { - /* only write if 3-phase solution */ - wellrates.set( rt::wat, wv[ wellrate_index + 0 ] ); - wellrates.set( rt::oil, wv[ wellrate_index + 1 ] ); - wellrates.set( rt::gas, wv[ wellrate_index + 2 ] ); + for( int i = 0; i < np; ++i ) { + comp.rates.set( phs[ i ], *(rates + i) ); + } } - - const double bhp = this->bhp()[ w ]; - const double thp = this->thp()[ w ]; - - res.emplace( wells_->name[ w ], - data::Well { wellrates, bhp, thp, std::move( completions ) } ); + assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]); } return res;