Output info numerical problem with multisegment wells.

opm/simulators/wells/MultisegmentWell.hpp

@@ -348,6 +348,9 @@ namespace Opm {
                        DeferredLogger& deferred_logger) const override;
 
         FSInfo getFirstPerforationFluidStateInfo(const Simulator& simulator) const;
+
+        void outputDebugInfoNumericalProblem(const SingleWellState<Scalar, IndexTraits>& ws,
+                                             DeferredLogger& deferred_logger) const;
     };
 
 }

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -574,6 +574,15 @@ getResidualMeasureValue(const WellState<Scalar, IndexTraits>& well_state,
     return sum;
 }
 
+template<typename FluidSystem, typename Indices>
+std::string
+MultisegmentWellEval<FluidSystem,Indices>::
+debugInfo() const
+{
+    // only outputing the primary variables for now
+    return primary_variables_.debugInfo();
+}
+
 #include <opm/simulators/utils/InstantiationIndicesMacros.hpp>
 
 INSTANTIATE_TYPE_INDICES(MultisegmentWellEval, double)

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -136,6 +136,8 @@ class MultisegmentWellEval : public MultisegmentWellGeneric<typename FluidSystem
 
     const WellInterfaceIndices<FluidSystem,Indices>& baseif_;
 
+    std::string debugInfo() const;
+
     Equations linSys_; //!< The equation system
     PrimaryVariables primary_variables_; //!< The primary variables
     MSWSegments segments_; //!< Segment properties

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -702,6 +702,27 @@ outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const
     }
 }
 
+template<class FluidSystem, class Indices>
+std::string
+MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
+debugInfo() const
+{
+    std::string msg = fmt::format("Well {} primary variables:\n", this->well_.name());
+    for (std::size_t seg = 0; seg < value_.size(); ++seg) {
+        fmt::format_to(std::back_inserter(msg), "  Segment {:4}: ", seg);
+        fmt::format_to(std::back_inserter(msg), "  SPres: {:8.2e} bar, ", value_[seg][SPres] / unit::barsa);
+        fmt::format_to(std::back_inserter(msg), "  WQTotal: {: 8.2e} m3/s, ", value_[seg][WQTotal]);
+        if (has_wfrac_variable) {
+            fmt::format_to(std::back_inserter(msg), "  WFrac: {:5.3f}, ", value_[seg][WFrac]);
+        }
+        if (has_gfrac_variable) {
+            fmt::format_to(std::back_inserter(msg), "  GFrac: {:5.3f},", value_[seg][GFrac]);
+        }
+        fmt::format_to(std::back_inserter(msg), "\n");
+    }
+    return msg;
+}
+
 #include <opm/simulators/utils/InstantiationIndicesMacros.hpp>
 
 INSTANTIATE_TYPE_INDICES(MultisegmentWellPrimaryVariables, double)

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -146,6 +146,9 @@ class MultisegmentWellPrimaryVariables
     //! output the segments with pressure close to lower pressure limit for debugging purpose
     void outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const;
 
+    //! output the debug information of the primary variables for debugging purpose
+    std::string debugInfo() const;
+
 private:
     //! \brief Initialize evaluations from values.
     void setEvaluationsFromValues();

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -283,6 +283,8 @@ namespace Opm
             // this is the process with rank zero)
             deferred_logger.problem("In MultisegmentWell::recoverWellSolutionAndUpdateWellState for well "
                                     + this->name() +": "+exp.what());
+            const auto& ws = well_state.well(this->index_of_well_);
+            this->outputDebugInfoNumericalProblem(ws, deferred_logger);
             throw;
         }
     }
@@ -632,6 +634,8 @@ namespace Opm
             // this is the process with rank zero)
             deferred_logger.problem("In MultisegmentWell::solveEqAndUpdateWellState for well "
                                     + this->name() +": "+exp.what());
+            const auto& ws = well_state.well(this->index_of_well_);
+            this->outputDebugInfoNumericalProblem(ws, deferred_logger);
             throw;
         }
     }
@@ -1607,6 +1611,8 @@ namespace Opm
                 // this is the process with rank zero)
                 deferred_logger.problem("In MultisegmentWell::iterateWellEqWithControl for well "
                                         + this->name() +": "+exp.what());
+                const auto& ws = well_state.well(this->index_of_well_);
+                this->outputDebugInfoNumericalProblem(ws, deferred_logger);
                 throw;
             }
         }
@@ -1788,6 +1794,8 @@ namespace Opm
                 // this is the process with rank zero)
                 deferred_logger.problem("In MultisegmentWell::iterateWellEqWithSwitching for well "
                                         + this->name() +": "+exp.what());
+                const auto& ws = well_state.well(this->index_of_well_);
+                this->outputDebugInfoNumericalProblem(ws, deferred_logger);
                 throw;
             }
         }
@@ -2379,6 +2387,19 @@ namespace Opm
         return this->parallel_well_info_.communication().size() == 1 ? info : this->pw_info_.broadcastFirstPerforationValue(info);
     }
 
+    template <typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    outputDebugInfoNumericalProblem(const SingleWellState<Scalar, IndexTraits>& ws,
+                                    DeferredLogger& deferred_logger) const
+    {
+        const std::string msg = ws.briefDebugInfo()
+                              + this->primary_variables_.debugInfo()
+                              + ws.segmentsDebugInfo()
+                              + ws.connectionDebugInfo();
+        deferred_logger.debug(msg);
+    }
+
 } // namespace Opm
 
 #endif

opm/simulators/wells/SegmentState.cpp

@@ -21,6 +21,8 @@
 #include "config.h"
 #endif // HAVE_CONFIG_H
 
+#include "opm/input/eclipse/Units/Units.hpp"
+
 #include <opm/simulators/wells/SegmentState.hpp>
 
 #include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
@@ -31,6 +33,8 @@
 #include <stdexcept>
 #include <vector>
 
+#include <fmt/format.h>
+
 namespace {
 
 std::vector<int> make_segment_number(const Opm::WellSegments& segments)
@@ -147,6 +151,27 @@ bool SegmentState<Scalar>::operator==(const SegmentState& rhs) const
            this->m_segment_number == rhs.m_segment_number;
 }
 
+template<class Scalar>
+std::string
+SegmentState<Scalar>::debugInfo() const
+{
+    if (this->empty()) {
+        return "";
+    }
+    const std::size_t num_phases = this->rates.size() / this->size();
+    std::string info = "SegmentState:\n";
+    for (std::size_t i = 0; i < this->size(); ++i) {
+        info += fmt::format("  Segment {:4}: Segment number {:4}, Pressure = {:8.2e} bar, Rate {m3/s}: ",
+                            i, this->m_segment_number[i], this->pressure[i] / unit::barsa);
+        for (std::size_t p = 0; p < num_phases; ++p) {
+            info += fmt::format(" {: 8.2e}", this->rates[i * num_phases + p]);
+        }
+        info += fmt::format("\n");
+    }
+    return info;
+}
+
+
 template class SegmentState<double>;
 
 #if FLOW_INSTANTIATE_FLOAT

opm/simulators/wells/SegmentState.hpp

@@ -22,6 +22,7 @@
 
 #include <cstddef>
 #include <vector>
+#include <string>
 
 namespace Opm {
 class WellSegments;
@@ -104,6 +105,7 @@ class SegmentState
     std::vector<Scalar> pressure_drop_hydrostatic;
     std::vector<Scalar> pressure_drop_accel;
 
+    std::string debugInfo() const;
 private:
     std::vector<int>    m_segment_number;
 };

opm/simulators/wells/SingleWellState.cpp

@@ -28,6 +28,8 @@
 
 #include <opm/simulators/wells/PerforationData.hpp>
 
+#include <fmt/format.h>
+
 namespace Opm {
 
 template<typename Scalar, typename IndexTraits>
@@ -415,6 +417,58 @@ bool SingleWellState<Scalar, IndexTraits>::operator==(const SingleWellState& rhs
            this->group_target == rhs.group_target;
 }
 
+template<typename Scalar, typename IndexTraits>
+std::string
+SingleWellState<Scalar, IndexTraits>::
+segmentsDebugInfo() const
+{
+    return this->segments.debugInfo();
+}
+
+template<typename Scalar, typename IndexTraits>
+std::string
+SingleWellState<Scalar, IndexTraits>::
+connectionDebugInfo() const
+{
+    std::string info {"  Connection pressures and connection rates:\n"};
+    for (std::size_t perf = 0; perf < this->perf_data.size(); perf++) {
+        fmt::format_to(std::back_inserter(info),
+                       "  Connection {:4}: Pressure: {:8.2e} bar, Rate:",
+                       perf,
+                       this->perf_data.pressure[perf] / unit::barsa);
+        const auto& connection_rates = this->perf_data.phase_rates;
+        const std::size_t num_phases = this->pu.numActivePhases();
+
+        for (std::size_t p = 0; p < num_phases; ++p) {
+            fmt::format_to(std::back_inserter(info), " {: 8.2e} m3/s", connection_rates[perf * num_phases + p]);
+        }
+        fmt::format_to(std::back_inserter(info), "\n");
+    }
+
+    return info;
+}
+
+template<typename Scalar, typename IndexTraits>
+std::string
+SingleWellState<Scalar, IndexTraits>::
+briefDebugInfo() const
+{
+    std::string info = "Well name: " + this->name + " well state:\n";
+    fmt::format_to(std::back_inserter(info), " type: {}, staus: {}, control type: {}, BHP: {:8.2e} Bar, THP: {:8.2e} Bar\n",
+                   this->producer? "Producer" : " Injector", WellStatus2String(this->status),
+                   this->producer ? WellProducerCMode2String(this->production_cmode)
+                                  : WellInjectorCMode2String(this->injection_cmode),
+                   this->bhp / unit::barsa,
+                   this->thp / unit::barsa);
+
+    fmt::format_to(std::back_inserter(info), "  Surface rates (m3/s): ");
+    for (const auto& r : this->surface_rates) {
+        fmt::format_to(std::back_inserter(info), " {:8.2e}", r);
+    }
+    fmt::format_to(std::back_inserter(info), "\n");
+    return info;
+}
+
 template class SingleWellState<double, BlackOilDefaultFluidSystemIndices>;
 
 #if FLOW_INSTANTIATE_FLOAT

opm/simulators/wells/SingleWellState.hpp

@@ -169,6 +169,10 @@ class SingleWellState {
     Scalar sum_filtrate_rate() const;
     Scalar sum_filtrate_total() const;
 
+    std::string briefDebugInfo() const;
+    std::string segmentsDebugInfo() const;
+    std::string connectionDebugInfo() const;
+
 private:
     Scalar sum_connection_rates(const std::vector<Scalar>& connection_rates) const;
 };