Merge pull request #1602 from atgeirr/convergence-status

Add convergence status class

CMakeLists_files.cmake

@@ -131,6 +131,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_blackoil_amg.cpp
   tests/test_block.cpp
   tests/test_boprops_ad.cpp
+  tests/test_convergencereport.cpp
   tests/test_graphcoloring.cpp
   tests/test_rateconverter.cpp
   tests/test_span.cpp
@@ -442,6 +443,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/timestepping/AdaptiveTimeStepping.hpp
   opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
   opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp
+  opm/simulators/timestepping/ConvergenceReport.hpp
   opm/simulators/timestepping/TimeStepControl.hpp
   opm/simulators/timestepping/TimeStepControlInterface.hpp
   opm/simulators/timestepping/SimulatorTimer.hpp

opm/autodiff/BlackoilWellModel.hpp

@@ -234,8 +234,6 @@ namespace Opm {
             // a vector of all the wells.
             std::vector<WellInterfacePtr > well_container_;
 
-            using ConvergenceReport = typename WellInterface<TypeTag>::ConvergenceReport;
-
             // create the well container
             std::vector<WellInterfacePtr > createWellContainer(const int time_step);
 

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -665,40 +665,48 @@ namespace Opm {
         for (const auto& well : well_container_) {
             report += well->getWellConvergence(B_avg);
         }
+        ConvergenceReport::Severity severity = report.severityOfWorstFailure();
 
         // checking NaN residuals
         {
-            bool nan_residual_found = report.nan_residual_found;
+            // Debug reporting.
+            for (const auto& f : report.wellFailures()) {
+                if (f.severity() == ConvergenceReport::Severity::NotANumber) {
+                    OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
+                }
+            }
+
+            // Throw if any nan residual found.
+            bool nan_residual_found = (severity == ConvergenceReport::Severity::NotANumber);
             const auto& grid = ebosSimulator_.vanguard().grid();
             int value = nan_residual_found ? 1 : 0;
-
             nan_residual_found = grid.comm().max(value);
-
             if (nan_residual_found) {
-                for (const auto& well : report.nan_residual_wells) {
-                    OpmLog::debug("NaN residual found with phase " + well.phase_name + " for well " + well.well_name);
-                }
                 OPM_THROW(Opm::NumericalIssue, "NaN residual found!");
             }
         }
 
         // checking too large residuals
         {
-            bool too_large_residual_found = report.too_large_residual_found;
+            // Debug reporting.
+            for (const auto& f : report.wellFailures()) {
+                if (f.severity() == ConvergenceReport::Severity::TooLarge) {
+                    OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
+                }
+            }
+
+            // Throw if any too large residual found.
+            bool too_large_residual_found = (severity == ConvergenceReport::Severity::TooLarge);
             const auto& grid = ebosSimulator_.vanguard().grid();
             int value = too_large_residual_found ? 1 : 0;
-
             too_large_residual_found = grid.comm().max(value);
             if (too_large_residual_found) {
-                for (const auto& well : report.too_large_residual_wells) {
-                    OpmLog::debug("Too large residual found with phase " + well.phase_name + " fow well " + well.well_name);
-                }
                 OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
             }
         }
 
         // checking convergence
-        bool converged_well = report.converged;
+        bool converged_well = report.converged();
         {
             const auto& grid = ebosSimulator_.vanguard().grid();
             int value = converged_well ? 1 : 0;

opm/autodiff/MultisegmentWell.hpp

@@ -71,7 +71,6 @@ namespace Opm
         static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq : numEq + 1;
 
         using typename Base::Scalar;
-        using typename Base::ConvergenceReport;
 
         /// the matrix and vector types for the reservoir
         using typename Base::Mat;

opm/autodiff/MultisegmentWell_impl.hpp

@@ -405,7 +405,7 @@ namespace Opm
 
 
     template <typename TypeTag>
-    typename MultisegmentWell<TypeTag>::ConvergenceReport
+    ConvergenceReport
     MultisegmentWell<TypeTag>::
     getWellConvergence(const std::vector<double>& B_avg) const
     {
@@ -419,62 +419,79 @@ namespace Opm
             }
         }
 
+        using CR = ConvergenceReport;
+        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
+        switch(well_controls_get_current_type(well_controls_)) {
+            case THP:
+                ctrltype = CR::WellFailure::Type::ControlTHP;
+                break;
+            case BHP:
+                ctrltype = CR::WellFailure::Type::ControlBHP;
+                break;
+            case RESERVOIR_RATE:
+            case SURFACE_RATE:
+                ctrltype = CR::WellFailure::Type::ControlRate;
+                break;
+            default:
+                OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name());
+        }
+        assert(ctrltype != CR::WellFailure::Type::Invalid);
+
         std::vector<double> maximum_residual(numWellEq, 0.0);
 
         ConvergenceReport report;
+        const int dummy_component = -1;
         // TODO: the following is a little complicated, maybe can be simplified in some way?
-        for (int seg = 0; seg < numberOfSegments(); ++seg) {
+        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
-            for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
+            for (int seg = 0; seg < numberOfSegments(); ++seg) {
                 if (eq_idx < num_components_) { // phase or component mass equations
                     const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
-                    // TODO: the report can not handle the segment number yet.
+                    if (flux_residual > maximum_residual[eq_idx]) {
-                    if (std::isnan(flux_residual)) {
+                        maximum_residual[eq_idx] = flux_residual;
-                        report.nan_residual_found = true;
-                        const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx));
-                        const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
-                        report.nan_residual_wells.push_back(problem_well);
-                    } else if (flux_residual > param_.max_residual_allowed_) {
-                        report.too_large_residual_found = true;
-                        const auto& compName = FluidSystem::componentName(Indices::activeToCanonicalComponentIndex(eq_idx));
-                        const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
-                        report.nan_residual_wells.push_back(problem_well);
-                    } else { // it is a normal residual
-                        if (flux_residual > maximum_residual[eq_idx]) {
-                            maximum_residual[eq_idx] = flux_residual;
-                        }
                     }
-                } else { // pressure equation
+                } else { // pressure or control equation
-                    // TODO: we should distinguish the rate control equations, bhp control equations
+                    if (seg == 0) {
-                    // and the oridnary pressure equations
+                        // Control equation
-                    const double pressure_residal = abs_residual[seg][eq_idx];
+                        const double control_residual = abs_residual[seg][eq_idx];
-                    const std::string eq_name("Pressure");
+                        if (std::isnan(control_residual)) {
-                    if (std::isnan(pressure_residal)) {
+                            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
-                        report.nan_residual_found = true;
+                        } else if (control_residual > param_.max_residual_allowed_) {
-                        const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
+                            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
-                        report.nan_residual_wells.push_back(problem_well);
+                        } else if (control_residual > param_.tolerance_wells_) {
-                    } else if (std::isinf(pressure_residal)) {
+                            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
-                        report.too_large_residual_found = true;
+                        }
-                        const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name};
+                    } else {
-                        report.nan_residual_wells.push_back(problem_well);
+                        // Pressure equation
-                    } else { // it is a normal residual
+                        const double pressure_residual = abs_residual[seg][eq_idx];
-                        if (pressure_residal > maximum_residual[eq_idx]) {
+                        if (pressure_residual > maximum_residual[eq_idx]) {
-                            maximum_residual[eq_idx] = pressure_residal;
+                            maximum_residual[eq_idx] = pressure_residual;
                         }
                     }
                 }
             }
         }
 
-        if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found
+        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
-            // check convergence for flux residuals
+            if (eq_idx < num_components_) { // phase or component mass equations
-            for ( int comp_idx = 0; comp_idx < num_components_; ++comp_idx)
+                const double flux_residual = maximum_residual[eq_idx];
-            {
+                // TODO: the report can not handle the segment number yet.
-                report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_);
+                if (std::isnan(flux_residual)) {
+                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, name()});
+                } else if (flux_residual > param_.max_residual_allowed_) {
+                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, name()});
+                } else if (flux_residual > param_.tolerance_wells_) {
+                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()});
+                }
+            } else { // pressure equation
+                const double pressure_residual = maximum_residual[eq_idx];
+                if (std::isnan(pressure_residual)) {
+                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()});
+                } else if (std::isinf(pressure_residual)) {
+                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, name()});
+                } else if (pressure_residual > param_.tolerance_pressure_ms_wells_) {
+                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()});
+                }
             }
-
-            report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_);
-        } else { // abnormal values found and no need to check the convergence
-            report.converged = false;
         }
 
         return report;
@@ -1726,9 +1743,8 @@ namespace Opm
             // const std::vector<double> B {0.8, 0.8, 0.008};
             const std::vector<double> B {0.5, 0.5, 0.005};
 
-            const ConvergenceReport report = getWellConvergence(B);
+            const auto report = getWellConvergence(B);
-
+            if (report.converged()) {
-            if (report.converged) {
                 break;
             }
 

opm/autodiff/StandardWell.hpp

@@ -86,7 +86,6 @@ namespace Opm
         static const int Bhp = numWellEq - numWellControlEq;
 
         using typename Base::Scalar;
-        using typename Base::ConvergenceReport;
 
 
         using Base::name;

opm/autodiff/StandardWell_impl.hpp

@@ -1483,7 +1483,7 @@ namespace Opm
 
 
     template<typename TypeTag>
-    typename StandardWell<TypeTag>::ConvergenceReport
+    ConvergenceReport
     StandardWell<TypeTag>::
     getWellConvergence(const std::vector<double>& B_avg) const
     {
@@ -1509,6 +1509,8 @@ namespace Opm
         }
 
         ConvergenceReport report;
+        using CR = ConvergenceReport;
+        CR::WellFailure::Type type = CR::WellFailure::Type::MassBalance;
         // checking if any NaN or too large residuals found
         for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
@@ -1517,62 +1519,46 @@ namespace Opm
 
             const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
             const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
-            const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
+            const int compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
 
             if (std::isnan(well_flux_residual[compIdx])) {
-                report.nan_residual_found = true;
+                report.setWellFailed({type, CR::Severity::NotANumber, compIdx, name()});
-                const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
+            } else if (well_flux_residual[compIdx] > maxResidualAllowed) {
-                report.nan_residual_wells.push_back(problem_well);
+                report.setWellFailed({type, CR::Severity::TooLarge, compIdx, name()});
-            } else {
+            } else if (well_flux_residual[compIdx] > tol_wells) {
-                if (well_flux_residual[compIdx] > maxResidualAllowed) {
+                report.setWellFailed({type, CR::Severity::Normal, compIdx, name()});
-                    report.too_large_residual_found = true;
-                    const typename ConvergenceReport::ProblemWell problem_well = {name(), compName};
-                    report.too_large_residual_wells.push_back(problem_well);
-                }
             }
         }
 
-
         // processing the residual of the well control equation
         const double well_control_residual = res[numWellEq - 1];
         // TODO: we should have better way to specify the control equation tolerance
         double control_tolerance = 0.;
         switch(well_controls_get_current_type(well_controls_)) {
             case THP:
+                type = CR::WellFailure::Type::ControlTHP;
+                control_tolerance = 1.e3; // 0.01 bar
+                break;
             case BHP:  // pressure type of control
+                type = CR::WellFailure::Type::ControlBHP;
                 control_tolerance = 1.e3; // 0.01 bar
                 break;
             case RESERVOIR_RATE:
             case SURFACE_RATE:
+                type = CR::WellFailure::Type::ControlRate;
                 control_tolerance = 1.e-4; // smaller tolerance for rate control
                 break;
             default:
                 OPM_THROW(std::runtime_error, "Unknown well control control types for well " << name());
         }
 
-        const bool control_eq_converged = well_control_residual < control_tolerance;
+        const int dummy_component = -1;
-
         if (std::isnan(well_control_residual)) {
-            report.nan_residual_found = true;
+            report.setWellFailed({type, CR::Severity::NotANumber, dummy_component, name()});
-            const typename ConvergenceReport::ProblemWell problem_well = {name(), "control"};
+        } else if (well_control_residual > maxResidualAllowed * 10.) {
-            report.nan_residual_wells.push_back(problem_well);
+            report.setWellFailed({type, CR::Severity::TooLarge, dummy_component, name()});
-        } else {
+        } else if ( well_control_residual > control_tolerance) {
-            // TODO: for pressure control equations, it can be pretty big during Newton iteration
+            report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
-            if (well_control_residual > maxResidualAllowed * 10.) {
-                report.too_large_residual_found = true;
-                const typename ConvergenceReport::ProblemWell problem_well = {name(), "control"};
-                report.too_large_residual_wells.push_back(problem_well);
-            }
-        }
-
-        if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found
-            // check convergence
-            for ( int compIdx = 0; compIdx < num_components_; ++compIdx )
-            {
-                report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells) && control_eq_converged;
-            }
-        } else { // abnormal values found and no need to check the convergence
-            report.converged = false;
         }
 
         return report;

opm/autodiff/WellInterface.hpp

@@ -44,6 +44,7 @@
 #include <opm/autodiff/BlackoilModelParametersEbos.hpp>
 #include <opm/autodiff/RateConverter.hpp>
 
+#include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
 #include<dune/common/fmatrix.hh>
@@ -137,37 +138,6 @@ namespace Opm
 
         virtual void initPrimaryVariablesEvaluation() const = 0;
 
-        /// a struct to collect information about the convergence checking
-        struct ConvergenceReport {
-            struct ProblemWell {
-                std::string well_name;
-                std::string phase_name;
-            };
-            bool converged = true;
-            bool nan_residual_found = false;
-            std::vector<ProblemWell> nan_residual_wells;
-            // We consider Inf is large residual here
-            bool too_large_residual_found = false;
-            std::vector<ProblemWell> too_large_residual_wells;
-
-            ConvergenceReport& operator+=(const ConvergenceReport& rhs) {
-                converged = converged && rhs.converged;
-                nan_residual_found = nan_residual_found || rhs.nan_residual_found;
-                if (rhs.nan_residual_found) {
-                    for (const ProblemWell& well : rhs.nan_residual_wells) {
-                        nan_residual_wells.push_back(well);
-                    }
-                }
-                too_large_residual_found = too_large_residual_found || rhs.too_large_residual_found;
-                if (rhs.too_large_residual_found) {
-                    for (const ProblemWell& well : rhs.too_large_residual_wells) {
-                        too_large_residual_wells.push_back(well);
-                    }
-                }
-                return *this;
-            }
-        };
-
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const = 0;
 
         virtual void solveEqAndUpdateWellState(WellState& well_state) = 0;

opm/autodiff/WellInterface_impl.hpp

@@ -932,10 +932,8 @@ namespace Opm
         do {
             assembleWellEq(ebosSimulator, dt, well_state, true);
 
-            ConvergenceReport report;
+            auto report = getWellConvergence(B_avg);
-            report = getWellConvergence(B_avg);
+            converged = report.converged();
-            converged = report.converged;
-
             if (converged) {
                 break;
             }

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -0,0 +1,176 @@
+/*
+  Copyright 2018 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2018 Equinor.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_CONVERGENCEREPORT_HEADER_INCLUDED
+#define OPM_CONVERGENCEREPORT_HEADER_INCLUDED
+
+#include <cassert>
+#include <numeric>
+#include <string>
+#include <vector>
+
+namespace Opm
+{
+
+    /// Represents the convergence status of the whole simulator, to
+    /// make it possible to query and store the reasons for
+    /// convergence failures.
+    class ConvergenceReport
+    {
+    public:
+
+        // ----------- Types -----------
+
+        enum Status { AllGood            = 0,
+                      ReservoirFailed    = 1 << 0,
+                      WellFailed         = 1 << 1 };
+        enum struct Severity { None       = 0,
+                               Normal     = 1,
+                               TooLarge   = 2,
+                               NotANumber = 3 };
+        class ReservoirFailure
+        {
+        public:
+            enum struct Type { Invalid, MassBalance, Cnv };
+            ReservoirFailure(Type t, Severity s, int phase, int cell_index)
+                : type_(t), severity_(s), phase_(phase), cell_index_(cell_index)
+            {
+            }
+            Type type() const { return type_; }
+            Severity severity() const { return severity_; }
+            int phase() const { return phase_; }
+            int cellIndex() const { return cell_index_; }
+        private:
+            Type type_;
+            Severity severity_;
+            int phase_;
+            int cell_index_;
+        };
+        class WellFailure
+        {
+        public:
+            enum struct Type { Invalid, MassBalance, Pressure, ControlBHP, ControlTHP, ControlRate };
+            WellFailure(Type t, Severity s, int phase, const std::string& well_name)
+                : type_(t), severity_(s), phase_(phase), well_name_(well_name)
+            {
+            }
+            Type type() const { return type_; }
+            Severity severity() const { return severity_; }
+            int phase() const { return phase_; }
+            const std::string& wellName() const { return well_name_; }
+        private:
+            Type type_;
+            Severity severity_;
+            int phase_;
+            std::string well_name_;
+        };
+
+        // ----------- Mutating member functions -----------
+
+        ConvergenceReport()
+            : status_{AllGood}
+            , res_failures_{}
+            , well_failures_{}
+        {
+        }
+
+        void clear()
+        {
+            status_ = AllGood;
+            res_failures_.clear();
+            well_failures_.clear();
+        }
+
+        void setReservoirFailed(const ReservoirFailure& rf)
+        {
+            status_ = static_cast<Status>(status_ | ReservoirFailed);
+            res_failures_.push_back(rf);
+        }
+
+        void setWellFailed(const WellFailure& wf)
+        {
+            status_ = static_cast<Status>(status_ | WellFailed);
+            well_failures_.push_back(wf);
+        }
+
+        ConvergenceReport& operator+=(const ConvergenceReport& other)
+        {
+            status_ = static_cast<Status>(status_ | other.status_);
+            res_failures_.insert(res_failures_.end(), other.res_failures_.begin(), other.res_failures_.end());
+            well_failures_.insert(well_failures_.end(), other.well_failures_.begin(), other.well_failures_.end());
+            assert(reservoirFailed() != res_failures_.empty());
+            assert(wellFailed() != well_failures_.empty());
+            return *this;
+        }
+
+        // ----------- Const member functions (queries) -----------
+
+        bool converged() const
+        {
+            return status_ == AllGood;
+        }
+
+        bool reservoirFailed() const
+        {
+            return status_ & ReservoirFailed;
+        }
+
+        bool wellFailed() const
+        {
+            return status_ & WellFailed;
+        }
+
+        const std::vector<ReservoirFailure>& reservoirFailures() const
+        {
+            return res_failures_;
+        }
+
+        const std::vector<WellFailure>& wellFailures() const
+        {
+            return well_failures_;
+        }
+
+        Severity severityOfWorstFailure() const
+        {
+            // A function to get the worst of two severities.
+            auto smax = [](Severity s1, Severity s2) {
+                return s1 < s2 ? s2 : s1;
+            };
+            auto s = Severity::None;
+            for (const auto f : res_failures_) {
+                s = smax(s, f.severity());
+            }
+            for (const auto f : well_failures_) {
+                s = smax(s, f.severity());
+            }
+            return s;
+        }
+
+    private:
+
+        // ----------- Member variables -----------
+        Status status_;
+        std::vector<ReservoirFailure> res_failures_;
+        std::vector<WellFailure> well_failures_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_CONVERGENCEREPORT_HEADER_INCLUDED

tests/test_convergencereport.cpp

@@ -0,0 +1,131 @@
+/*
+  Copyright 2018 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2018 Equinor.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#define BOOST_TEST_MODULE ConvergenceReportTest
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/timestepping/ConvergenceReport.hpp>
+
+using CR = Opm::ConvergenceReport;
+
+BOOST_AUTO_TEST_CASE(DefaultConstructor)
+{
+    Opm::ConvergenceReport s;
+    BOOST_CHECK(s.converged());
+    BOOST_CHECK(!s.reservoirFailed());
+    BOOST_CHECK(!s.wellFailed());
+    BOOST_CHECK(s.severityOfWorstFailure() == CR::Severity::None);
+}
+
+BOOST_AUTO_TEST_CASE(Failures)
+{
+    Opm::ConvergenceReport s1;
+    s1.setReservoirFailed({CR::ReservoirFailure::Type::Cnv, CR::Severity::Normal, 2, 100});
+    {
+        BOOST_CHECK(!s1.converged());
+        BOOST_CHECK(s1.reservoirFailed());
+        BOOST_CHECK(!s1.wellFailed());
+        BOOST_REQUIRE(s1.reservoirFailures().size() == 1);
+        const auto f = s1.reservoirFailures()[0];
+        BOOST_CHECK(f.type() == CR::ReservoirFailure::Type::Cnv);
+        BOOST_CHECK(f.severity() == CR::Severity::Normal);
+        BOOST_CHECK(f.phase() == 2);
+        BOOST_CHECK(f.cellIndex() == 100);
+        BOOST_CHECK(s1.wellFailures().empty());
+        BOOST_CHECK(s1.severityOfWorstFailure() == CR::Severity::Normal);
+    }
+
+    Opm::ConvergenceReport s2;
+    s2.setWellFailed({CR::WellFailure::Type::ControlTHP, CR::Severity::Normal, -1, "PRODUCER-123"});
+    s2.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, 2, "INJECTOR-XYZ"});
+    {
+        BOOST_CHECK(!s2.converged());
+        BOOST_CHECK(!s2.reservoirFailed());
+        BOOST_CHECK(s2.wellFailed());
+        BOOST_CHECK(s2.reservoirFailures().empty());
+        BOOST_REQUIRE(s2.wellFailures().size() == 2);
+        const auto f0 = s2.wellFailures()[0];
+        BOOST_CHECK(f0.type() == CR::WellFailure::Type::ControlTHP);
+        BOOST_CHECK(f0.severity() == CR::Severity::Normal);
+        BOOST_CHECK(f0.phase() == -1);
+        BOOST_CHECK(f0.wellName() == "PRODUCER-123");
+        const auto f1 = s2.wellFailures()[1];
+        BOOST_CHECK(f1.type() == CR::WellFailure::Type::MassBalance);
+        BOOST_CHECK(f1.severity() == CR::Severity::TooLarge);
+        BOOST_CHECK(f1.phase() == 2);
+        BOOST_CHECK(f1.wellName() == "INJECTOR-XYZ");
+        BOOST_CHECK(s2.severityOfWorstFailure() == CR::Severity::TooLarge);
+    }
+
+    s1 += s2;
+    {
+        BOOST_CHECK(!s1.converged());
+        BOOST_CHECK(s1.reservoirFailed());
+        BOOST_CHECK(s1.wellFailed());
+        BOOST_REQUIRE(s1.reservoirFailures().size() == 1);
+        const auto f = s1.reservoirFailures()[0];
+        BOOST_CHECK(f.type() == CR::ReservoirFailure::Type::Cnv);
+        BOOST_CHECK(f.severity() == CR::Severity::Normal);
+        BOOST_CHECK(f.phase() == 2);
+        BOOST_CHECK(f.cellIndex() == 100);
+        BOOST_REQUIRE(s1.wellFailures().size() == 2);
+        const auto f0 = s1.wellFailures()[0];
+        BOOST_CHECK(f0.type() == CR::WellFailure::Type::ControlTHP);
+        BOOST_CHECK(f0.severity() == CR::Severity::Normal);
+        BOOST_CHECK(f0.phase() == -1);
+        BOOST_CHECK(f0.wellName() == "PRODUCER-123");
+        const auto f1 = s1.wellFailures()[1];
+        BOOST_CHECK(f1.type() == CR::WellFailure::Type::MassBalance);
+        BOOST_CHECK(f1.severity() == CR::Severity::TooLarge);
+        BOOST_CHECK(f1.phase() == 2);
+        BOOST_CHECK(f1.wellName() == "INJECTOR-XYZ");
+        BOOST_CHECK(s1.severityOfWorstFailure() == CR::Severity::TooLarge);
+    }
+
+    s1.clear();
+    {
+        BOOST_CHECK(s1.converged());
+        BOOST_CHECK(!s1.reservoirFailed());
+        BOOST_CHECK(!s1.wellFailed());
+        BOOST_CHECK(s1.severityOfWorstFailure() == CR::Severity::None);
+    }
+
+    s1 += s2;
+    {
+        BOOST_CHECK(!s1.converged());
+        BOOST_CHECK(!s1.reservoirFailed());
+        BOOST_CHECK(s1.wellFailed());
+        BOOST_CHECK(s1.reservoirFailures().empty());
+        BOOST_REQUIRE(s1.wellFailures().size() == 2);
+        const auto f0 = s1.wellFailures()[0];
+        BOOST_CHECK(f0.type() == CR::WellFailure::Type::ControlTHP);
+        BOOST_CHECK(f0.severity() == CR::Severity::Normal);
+        BOOST_CHECK(f0.phase() == -1);
+        BOOST_CHECK(f0.wellName() == "PRODUCER-123");
+        const auto f1 = s1.wellFailures()[1];
+        BOOST_CHECK(f1.type() == CR::WellFailure::Type::MassBalance);
+        BOOST_CHECK(f1.severity() == CR::Severity::TooLarge);
+        BOOST_CHECK(f1.phase() == 2);
+        BOOST_CHECK(f1.wellName() == "INJECTOR-XYZ");
+        BOOST_CHECK(s1.severityOfWorstFailure() == CR::Severity::TooLarge);
+    }
+}
+