Merge remote-tracking branch 'origin/master' into frankenstein_merge_master_v7

* origin/master:
  Bugfix: use correct object in sequential model.
  Bugfix: properly handle two-phase cases.
  Bugfix: pass both normal and extra data to output.
  Update well data output integration.
  Adapt to API change in opm-parser.
  Ensure logging only on first rank.
  Add isIORank_ member to avoid repeated calls.
  Rename isRankZero() -> isIORank() for consistency.
  Ensure only first-rank logging.
  Add isRankZero() utility.
  Bugfix: missing return.
  do not set limits for prt log and set correct value for Note.
  correct function call order and add whitespace.
  supprot MESSAGES default vaule.
  wellToState reads new opm-output data exchange
  Write control to data::Wells. Missing ability to restore.
  Avoid using buggy wellstate api
  Stop report early if there are no wells
  Update to interface change in opm-output

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -1744,16 +1744,28 @@ typedef Eigen::Array<double,
             if (std::isnan(mass_balance_residual[idx])
                 || std::isnan(CNV[idx])
                 || (idx < np && std::isnan(well_flux_residual[idx]))) {
-                OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
+                const auto msg = std::string("NaN residual for phase ") +  materialName(idx);
+                if (terminal_output_) {
+                    OpmLog::problem(msg);
+                }
+                OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
             }
             if (mass_balance_residual[idx] > 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<double,
         // if one of the residuals is NaN, throw exception, so that the solver can be restarted
         for (int idx = 0; idx < np; ++idx) {
             if (std::isnan(well_flux_residual[idx])) {
-                OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
+                const auto msg = std::string("NaN residual for phase ") +  materialName(idx);
+                if (terminal_output_) {
+                    OpmLog::problem(msg);
+                }
+                OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
             }
             if (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);
             }
         }
 
@@ -2198,10 +2218,12 @@ typedef Eigen::Array<double,
         {
             //Accumulate phases for each region
             for (int phase = 0; phase < maxnp; ++phase) {
-                for (int c = 0; c < nc; ++c) {
+                if (active_[ phase ]) {
-                    const int region = fipnum[c] - 1;
+                    for (int c = 0; c < nc; ++c) {
-                    if (region != -1) {
+                        const int region = fipnum[c] - 1;
-                        values[region][phase] += sd_.fip[phase][c];
+                        if (region != -1) {
+                            values[region][phase] += sd_.fip[phase][c];
+                        }
                     }
                 }
             }

opm/autodiff/BlackoilSequentialModel.hpp

@@ -247,7 +247,7 @@ namespace Opm {
         /// Return the well model
         const WellModel& wellModel() const
         {
-            return pressure_model_->wellModel();
+            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();
         }
 
 

opm/autodiff/Compat.hpp

@@ -21,10 +21,13 @@
 #ifndef OPM_SIMULATORS_COMPAT_HPP
 #define OPM_SIMULATORS_COMPAT_HPP
 
+#include <algorithm>
+#include <cassert>
+
 #include <opm/common/data/SimulationDataContainer.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilSolventState.hpp>
 #include <opm/output/data/Cells.hpp>
 #include <opm/output/data/Solution.hpp>
@@ -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 ) {
+    const auto np = phases.num_phases;
-    state.bhp() = wells.bhp;
+
-    state.temperature() = wells.temperature;
+    std::vector< rt > phs( np );
-    state.wellRates() = wells.well_rate;
+    if( phases.phase_used[BlackoilPhases::Aqua] ) {
-    state.perfPress() = wells.perf_pressure;
+        phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
-    state.perfRates() = wells.perf_rate;
+    }
+
+    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;
+        }
+    }
 }
 
 }

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<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
             OpmLog::addBackend( "DEBUGLOG" ,  debugLog);
+            const auto& msgLimits = eclipse_state_->getSchedule().getMessageLimits();
+            const std::map<int64_t, int> 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<MessageLimiter>());
             prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
-            streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10));
+            streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
             streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(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);

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.

opm/autodiff/NewtonIterationUtilities.cpp

@@ -22,6 +22,7 @@
 
 #include <opm/autodiff/NewtonIterationUtilities.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
 #include <opm/common/ErrorMacros.hpp>
 
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
@@ -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<const ParallelISTLInformation&>(parallel_info);
+            return info.communicator().rank() == 0;
+        } else {
+            return true;
+        }
+#else
+        static_cast<void>(parallel_info); // Suppress unused argument warning.
+        return true;
+#endif
+    }
+
+
 } // namespace Opm
 

opm/autodiff/NewtonIterationUtilities.hpp

@@ -21,6 +21,7 @@
 #define OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <boost/any.hpp>
 #include <vector>
 
 namespace Opm
@@ -58,6 +59,8 @@ namespace Opm
                             Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
                             AutoDiffBlock<double>::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
 

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

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),

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
+            if( nw == 0 ) return res;
-            // exception.
+            const int np = pu.num_phases;
-            const int np = nw ? this->numPhases() : -1;
+
+
+            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 ) {
+            for( const auto& wt : this->wellMap() ) {
-                using rt = data::Rates::opt;
+                const auto w = wt.second[ 0 ];
-                std::map< size_t, data::Completion > completions;
+                auto& well = res.at( wt.first );
+                well.control = this->currentControls()[ w ];
 
-                // completions aren't supported yet
+                int local_comp_index = 0;
-                //const auto* begin = wells_->well_connpos + w;
+                for( auto& comp : well.completions ) {
-                //const auto* end = wells_->well_connpos + w + 1;
+                    const auto rates = this->perfPhaseRates().begin()
-                //for( auto* i = begin; i != end; ++i ) {
+                                     + (np * wt.second[ 1 ])
-                //    const auto perfrate = this->perfPhaseRates().begin() + *i;
+                                     + (np * local_comp_index);
-                //    data::Rates perfrates;
+                    ++local_comp_index;
-                //    perfrates.set( rt::wat, *(perfrate + 0) );
-                //    perfrates.set( rt::oil, *(perfrate + 1) );
-                //    perfrates.set( rt::gas, *(perfrate + 2) );
 
-                //    const size_t active_index = wells_->well_cells[ *i ];
+                    for( int i = 0; i < np; ++i ) {
-
+                        comp.rates.set( phs[ i ], *(rates + 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 ] );
                 }
-
+                assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
-                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 ) } );
             }
 
             return res;