fix some serious screw-ups

almost all of them were caused by recent changes in the master
branch:

- there were methods added which depend on the types `V` and
`DataBlock`. these do not make much sense in the context of the
frankenstein simulator. Also, these types are defined globally for the
whole Opm namespace in `BlackoilModelBase_impl.hpp` (which should be
prosecuted as a fellony IMO)! Besides this, their names are useless;
'V' is the letter which comes after `U` in the alphabet and when it
comes to computers basically everything can be seen as a chunk of data
(i.e., a `DataBlock`).
- it seems like the new and shiny dense-AD based well model was never
compiled with assertations enabled, at least some asserts referenced
non-existing variables.
- the recent output-related API changes were pretty unfortunate
because they had the effect of tying the (sub-optimal, IMO) internal
structure of the model even closer to the output code: as far as I can
see, `rq` does only make sense if the model works *exactly* like
BlackoilModelBase and friends. (for flow_ebos, this could be
replicated, but first it would be another unnecessary conversion step
and second, most of the quantities in `rq` are of type `ADB` and much
of the "frankenstein" excercise is devoted to getting rid of these.) I
thus reverted back to an old version of the output code and created a
`frankenstein` branch in my personal `opm-output` github fork.

CMakeLists_files.cmake

@@ -36,6 +36,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/autodiff/ImpesTPFAAD.cpp
   opm/autodiff/moduleVersion.cpp
   opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
+  opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp
   opm/autodiff/SimulatorIncompTwophaseAd.cpp
   opm/autodiff/TransportSolverTwophaseAd.cpp
   opm/autodiff/BlackoilPropsAdFromDeck.cpp

opm/autodiff/BlackoilModelBase.hpp

@@ -27,13 +27,13 @@
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
-#include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
-#include <opm/autodiff/IterationReport.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/autodiff/BlackoilModelEnums.hpp>
 #include <opm/autodiff/VFPProperties.hpp>
+#include <opm/autodiff/IterationReport.hpp>
+#include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
 

opm/autodiff/BlackoilModelEbos.hpp

@@ -42,6 +42,7 @@
 #include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
 #include <opm/autodiff/BlackoilDetails.hpp>
 #include <opm/autodiff/BlackoilModelEnums.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
@@ -52,6 +53,8 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/well_controls.h>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
@@ -221,7 +224,6 @@ namespace Opm {
                                            ReservoirState& reservoir_state,
                                            WellState& well_state)
         {
-            const double dt = timer.currentStepLength();
             if (iteration == 0) {
                 // For each iteration we store in a vector the norms of the residual of
                 // the mass balance for each active phase, the well flux and the well equations.
@@ -302,13 +304,6 @@ namespace Opm {
                 OPM_THROW(Opm::NumericalProblem,"no convergence");
             }
 
-            typedef double Scalar;
-            typedef Dune::FieldVector<Scalar, 3    >       VectorBlockType;
-            typedef Dune::FieldMatrix<Scalar, 3, 3 >      MatrixBlockType;
-            typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
-            typedef Dune::BlockVector<VectorBlockType>      BVector;
-
-            typedef Dune::MatrixAdapter<Mat,BVector,BVector> Operator;
             auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
 
@@ -927,6 +922,16 @@ namespace Opm {
             }
         }
 
+        std::vector<std::vector<double> >
+        computeFluidInPlace(const ReservoirState& x,
+                            const std::vector<int>& fipnum) const
+        {
+            OPM_THROW(std::logic_error,
+                      "computeFluidInPlace() not implemented by BlackoilModelEbos!");
+        }
+
+        const Simulator& ebosSimulator() const
+        { return ebosSimulator_; }
 
     protected:
 
@@ -1244,8 +1249,7 @@ namespace Opm {
         {
             const int numPhases = wells().number_of_phases;
             const int numCells = ebosJac.N();
-            const int cols = ebosJac.M();
+            assert( numCells == static_cast<int>(ebosJac.M()) );
-            assert( numCells == cols );
 
             // write the right-hand-side values from the ebosJac into the objects
             // allocated above.

opm/autodiff/BlackoilPressureModel.hpp

@@ -26,7 +26,6 @@
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
-#include <opm/autodiff/StandardWells.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
 
 #include <algorithm>

opm/autodiff/FlowMain.hpp

@@ -173,6 +173,7 @@ namespace Opm
         typedef BlackoilPropsAdFromDeck FluidProps;
         typedef FluidProps::MaterialLawManager MaterialLawManager;
         typedef typename Simulator::ReservoirState ReservoirState;
+        typedef typename Simulator::OutputWriter OutputWriter;
 
 
         // ------------   Data members   ------------
@@ -207,7 +208,7 @@ namespace Opm
         // distributeData()
         boost::any parallel_information_;
         // setupOutputWriter()
-        std::unique_ptr<BlackoilOutputWriter> output_writer_;
+        std::unique_ptr<OutputWriter> output_writer_;
         // setupLinearSolver
         std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
         // createSimulator()
@@ -673,11 +674,11 @@ namespace Opm
             // create output writer after grid is distributed, otherwise the parallel output
             // won't work correctly since we need to create a mapping from the distributed to
             // the global view
-            output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(),
+            output_writer_.reset(new OutputWriter(grid_init_->grid(),
-                                                          param_,
+                                                  param_,
-                                                          eclipse_state_,
+                                                  eclipse_state_,
-                                                          Opm::phaseUsageFromDeck(deck_),
+                                                  Opm::phaseUsageFromDeck(deck_),
-                                                          fluidprops_->permeability()));
+                                                  fluidprops_->permeability()));
         }
 
 

opm/autodiff/NewtonIterationBlackoilInterleaved.cpp

@@ -390,7 +390,7 @@ namespace Opm
             //const int np = residual.material_balance_eq.size();
             assert( np == int(residual.material_balance_eq.size()) );
             std::vector<ADB> eqs;
-            eqs.reserve(np + 1);
+            eqs.reserve(np + 2);
             for (int phase = 0; phase < np; ++phase) {
                 eqs.push_back(residual.material_balance_eq[phase]);
             }
@@ -401,14 +401,14 @@ namespace Opm
             if( hasWells )
             {
                 eqs.push_back(residual.well_flux_eq);
-                //eqs.push_back(residual.well_eq);
+                eqs.push_back(residual.well_eq);
 
                 // Eliminate the well-related unknowns, and corresponding equations.
-                elim_eqs.reserve(1);
+                elim_eqs.reserve(2);
                 elim_eqs.push_back(eqs[np]);
                 eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns.
-                //elim_eqs.push_back(eqs[np]);
+                elim_eqs.push_back(eqs[np]);
-                //eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
+                eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
                 assert(int(eqs.size()) == np);
             }
 
@@ -500,7 +500,7 @@ namespace Opm
             if ( hasWells ) {
                 // Compute full solution using the eliminated equations.
                 // Recovery in inverse order of elimination.
-                //dx = recoverVariable(elim_eqs[1], dx, np);
+                dx = recoverVariable(elim_eqs[1], dx, np);
                 dx = recoverVariable(elim_eqs[0], dx, np);
             }
             return dx;
@@ -575,7 +575,6 @@ namespace Opm
         computePressureIncrement(const LinearisedBlackoilResidual& residual)
         {
             typedef LinearisedBlackoilResidual::ADB ADB;
-            typedef ADB::V V;
 
             // Build the vector of equations (should be just a single material balance equation
             // in which the pressure equation is stored).

opm/autodiff/NewtonIterationUtilities.cpp

@@ -85,16 +85,9 @@ namespace Opm
         for (int eq = 0; eq < num_eq; ++eq) {
             jacs[eq].reserve(num_eq - 1);
             const std::vector<M>& Je = eqs[eq].derivative();
-            Sp Bb;
             const M& B = Je[n];
-            B.toSparse(Bb);
-            //std::cout << "B eigen" << std::endl;
-            //std::cout << Bb << std::endl;
-
             // Update right hand side.
             vals[eq] = eqs[eq].value().matrix() - B * Dibn;
-            //std::cout << "vals " << eq << std::endl;
-            //std::cout << vals[eq][0] << std::endl;
         }
         for (int var = 0; var < num_eq; ++var) {
             if (var == n) {
@@ -156,11 +149,6 @@ namespace Opm
         V equation_value = equation.value();
         ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs)));
         const M& C = eq_coll.derivative()[0];
-        typedef Eigen::SparseMatrix<double> Sp;
-        Sp Cc;
-        C.toSparse(Cc);
-        //std::cout << "C eigen" << std::endl;
-        //std::cout << Cc <<std::endl;
 
         // Use sparse LU to solve the block submatrices
         typedef Eigen::SparseMatrix<double> Sp;

opm/autodiff/NonlinearSolver.hpp

@@ -131,7 +131,11 @@ namespace Opm {
         /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
         std::vector<V>
         computeFluidInPlace(const ReservoirState& x,
-                            const std::vector<int>& fipnum) const;
+                            const std::vector<int>& fipnum) const
+        {
+            return model_->computeFluidInPlace(x, fipnum);
+        }
+
 
         /// Reference to physical model.
         const PhysicalModel& model() const;

opm/autodiff/NonlinearSolver_impl.hpp

@@ -24,6 +24,8 @@
 #define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
 
 #include <opm/autodiff/NonlinearSolver.hpp>
+#include <opm/common/Exceptions.hpp>
+#include <opm/common/ErrorMacros.hpp>
 
 namespace Opm
 {
@@ -99,15 +101,6 @@ namespace Opm
         return wellIterationsLast_;
     }
 
-    template <class PhysicalModel>
-    std::vector<V> 
-    NonlinearSolver<PhysicalModel>::computeFluidInPlace(const ReservoirState& x,
-                                                        const std::vector<int>& fipnum) const
-    {
-        return model_->computeFluidInPlace(x, fipnum);
-    }
-
-
     template <class PhysicalModel>
     int
     NonlinearSolver<PhysicalModel>::
@@ -119,6 +112,7 @@ namespace Opm
     }
 
 
+
     template <class PhysicalModel>
     int
     NonlinearSolver<PhysicalModel>::
@@ -177,6 +171,7 @@ namespace Opm
     }
 
 
+
     template <class PhysicalModel>
     void NonlinearSolver<PhysicalModel>::SolverParameters::
     reset()

opm/autodiff/SimFIBODetails.hpp

@@ -0,0 +1,151 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014-2016 IRIS AS
+  Copyright 2015 Andreas Lauser
+
+  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_SIM_FIBO_DETAILS_HPP
+#define OPM_SIM_FIBO_DETAILS_HPP
+
+#include <utility>
+#include <algorithm>
+#include <locale>
+#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
+#include <opm/core/utility/initHydroCarbonState.hpp>
+#include <opm/core/well_controls.h>
+#include <opm/core/wells/DynamicListEconLimited.hpp>
+
+namespace Opm
+{
+    namespace SimFIBODetails {
+        typedef std::unordered_map<std::string, const Well* > WellMap;
+
+        inline WellMap
+        mapWells(const std::vector< const Well* >& wells)
+        {
+            WellMap wmap;
+
+            for (std::vector< const Well* >::const_iterator
+                     w = wells.begin(), e = wells.end();
+                 w != e; ++w)
+            {
+                wmap.insert(std::make_pair((*w)->name(), *w));
+            }
+
+            return wmap;
+        }
+
+        inline int
+        resv_control(const WellControls* ctrl)
+        {
+            int i, n = well_controls_get_num(ctrl);
+
+            bool match = false;
+            for (i = 0; (! match) && (i < n); ++i) {
+                match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
+            }
+
+            if (! match) { i = 0; }
+
+            return i - 1; // -1 if no match, undo final "++" otherwise
+        }
+
+        inline bool
+        is_resv(const Wells& wells,
+                const int    w)
+        {
+            return (0 <= resv_control(wells.ctrls[w]));
+        }
+
+        inline bool
+        is_resv(const WellMap&     wmap,
+                const std::string& name,
+                const std::size_t  step)
+        {
+            bool match = false;
+
+            WellMap::const_iterator i = wmap.find(name);
+
+            if (i != wmap.end()) {
+                const Well* wp = i->second;
+
+                match = (wp->isProducer(step) &&
+                         wp->getProductionProperties(step)
+                         .hasProductionControl(WellProducer::RESV))
+                    ||  (wp->isInjector(step) &&
+                         wp->getInjectionProperties(step)
+                         .hasInjectionControl(WellInjector::RESV));
+            }
+
+            return match;
+        }
+
+        inline std::vector<int>
+        resvWells(const Wells*      wells,
+                  const std::size_t step,
+                  const WellMap&    wmap)
+        {
+            std::vector<int> resv_wells;
+            if( wells )
+            {
+                for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
+                    if (is_resv(*wells, w) ||
+                        ((wells->name[w] != 0) &&
+                         is_resv(wmap, wells->name[w], step)))
+                    {
+                        resv_wells.push_back(w);
+                    }
+                }
+            }
+
+            return resv_wells;
+        }
+
+        inline void
+        historyRates(const PhaseUsage&               pu,
+                     const WellProductionProperties& p,
+                     std::vector<double>&            rates)
+        {
+            assert (! p.predictionMode);
+            assert (rates.size() ==
+                    std::vector<double>::size_type(pu.num_phases));
+
+            if (pu.phase_used[ BlackoilPhases::Aqua ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Aqua ];
+
+                rates[i] = p.WaterRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Liquid ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Liquid ];
+
+                rates[i] = p.OilRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Vapour ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Vapour ];
+
+                rates[i] = p.GasRate;
+            }
+        }
+    } // namespace SimFIBODetails
+} // namespace Opm
+
+#endif // OPM_SIM_FIBO_DETAILS_HPP

opm/autodiff/SimulatorBase.hpp

@@ -21,14 +21,12 @@
 #ifndef OPM_SIMULATORBASE_HEADER_INCLUDED
 #define OPM_SIMULATORBASE_HEADER_INCLUDED
 
-#include <ewoms/common/timer.hh>
 #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
-#include <opm/autodiff/IterationReport.hpp>
-#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/common/ErrorMacros.hpp>
 
 #include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/BlackoilModel.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/RateConverter.hpp>

opm/autodiff/SimulatorBase_impl.hpp

@@ -26,6 +26,7 @@
 #include <opm/core/utility/initHydroCarbonState.hpp>
 #include <opm/core/well_controls.h>
 #include <opm/core/wells/DynamicListEconLimited.hpp>
+#include <opm/autodiff/BlackoilModel.hpp>
 
 namespace Opm
 {

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -18,17 +18,26 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
-#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
 
-#include <opm/autodiff/SimulatorBase.hpp>
+//#include <opm/autodiff/SimulatorBase.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp>
+#include <opm/autodiff/IterationReport.hpp>
 #include <opm/autodiff/NonlinearSolver.hpp>
 #include <opm/autodiff/BlackoilModelEbos.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/StandardWellsDense.hpp>
+#include <opm/autodiff/RateConverter.hpp>
+#include <opm/autodiff/SimFIBODetails.hpp>
 
+#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
+#include <opm/core/utility/initHydroCarbonState.hpp>
+#include <opm/core/utility/StopWatch.hpp>
 
+#include <opm/common/Exceptions.hpp>
+#include <opm/common/ErrorMacros.hpp>
 
 namespace Opm {
 
@@ -47,7 +56,7 @@ public:
 
     typedef WellStateFullyImplicitBlackoil WellState;
     typedef BlackoilState ReservoirState;
-    typedef BlackoilOutputWriter OutputWriter;
+    typedef BlackoilOutputWriterEbos OutputWriter;
     typedef BlackoilModelEbos Model;
     typedef BlackoilModelParameters ModelParameters;
     typedef NonlinearSolver<Model> Solver;
@@ -90,7 +99,7 @@ public:
                                        const bool has_disgas,
                                        const bool has_vapoil,
                                        std::shared_ptr<EclipseState> eclipse_state,
-                                       BlackoilOutputWriter& output_writer,
+                                       BlackoilOutputWriterEbos& output_writer,
                                        const std::vector<double>& threshold_pressures_by_face)
         : ebosSimulator_(ebosSimulator),
           param_(param),
@@ -166,15 +175,11 @@ public:
             adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
         }
 
-
-
-        output_writer_.writeInit( geo_.simProps(grid()) , geo_.nonCartesianConnections( ) );
-
         std::string restorefilename = param_.getDefault("restorefile", std::string("") );
         if( ! restorefilename.empty() )
         {
             // -1 means that we'll take the last report step that was written
-            const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
+//            const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
 
 //            output_writer_.restore( timer,
 //                                    state,
@@ -221,11 +226,6 @@ public:
             // give the polymer and surfactant simulators the chance to do their stuff
             handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
 
-            // write the inital state at the report stage
-            if (timer.initialStep()) {
-                output_writer_.writeTimeStep( timer, state, well_state );
-            }
-
             // Compute reservoir volumes for RESV controls.
             computeRESV(timer.currentStepNum(), wells, state, well_state);
 
@@ -237,6 +237,11 @@ public:
 
             auto solver = createSolver(well_model);
 
+            // write the inital state at the report stage
+            if (timer.initialStep()) {
+                output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
+            }
+
             if( terminal_output_ )
             {
                 std::ostringstream step_msg;
@@ -326,7 +331,7 @@ public:
             ++timer;
 
             // write simulation state at the report stage
-            output_writer_.writeTimeStep( timer, state, well_state );
+            output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
 
             prev_well_state = well_state;
             // The well potentials are only computed if they are needed

opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp

@@ -0,0 +1,282 @@
+/*
+  Copyright (c) 2014 SINTEF ICT, Applied Mathematics.
+  Copyright (c) 2015-2016 IRIS AS
+
+  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"
+
+#include "SimulatorFullyImplicitBlackoilOutputEbos.hpp"
+
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+#include <opm/output/Cells.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/core/utility/Compat.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/Units.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/BackupRestore.hpp>
+
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+//For OutputWriterHelper
+#include <map>
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+
+#ifdef HAVE_OPM_GRID
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <dune/common/version.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+#endif
+namespace Opm
+{
+    namespace detail {
+
+        struct WriterCallEbos : public ThreadHandle :: ObjectInterface
+        {
+            BlackoilOutputWriterEbos& writer_;
+            std::unique_ptr< SimulatorTimerInterface > timer_;
+            const SimulationDataContainer state_;
+            const WellState wellState_;
+            std::vector<data::CellData> simProps_;
+            const bool substep_;
+
+            explicit WriterCallEbos( BlackoilOutputWriterEbos& writer,
+                                 const SimulatorTimerInterface& timer,
+                                 const SimulationDataContainer& state,
+                                 const WellState& wellState,
+                                 const std::vector<data::CellData>& simProps,
+                                 bool substep )
+                : writer_( writer ),
+                  timer_( timer.clone() ),
+                  state_( state ),
+                  wellState_( wellState ),
+                  simProps_( simProps ),
+                  substep_( substep )
+            {
+            }
+
+            // callback to writer's serial writeTimeStep method
+            void run ()
+            {
+                // write data
+                writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ );
+            }
+        };
+    }
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepWithoutCellProperties(
+                  const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  bool substep)
+    {
+        std::vector<data::CellData> noCellProperties;
+        writeTimeStepWithCellProperties(timer, localState, localWellState, noCellProperties, substep);
+    }
+
+
+
+
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepWithCellProperties(
+                  const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  const std::vector<data::CellData>& cellData,
+                  bool substep)
+    {
+        bool isIORank = output_ ;
+        if( parallelOutput_ && parallelOutput_->isParallel() )
+        {
+            // If this is not the initial write and no substep, then the well
+            // state used in the computation is actually the one of the last
+            // step. We need that well state for the gathering. Otherwise
+            // It an exception with a message like "global state does not
+            // contain well ..." might be thrown.
+            int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
+                (timer.reportStepNum() - 1) : timer.reportStepNum();
+            // collect all solutions to I/O rank
+            isIORank = parallelOutput_->collectToIORank( localState, localWellState, wellStateStepNumber );
+        }
+
+        const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
+        const WellState& wellState  = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
+
+        // serial output is only done on I/O rank
+        if( isIORank )
+        {
+            if( asyncOutput_ ) {
+                // dispatch the write call to the extra thread
+                asyncOutput_->dispatch( detail::WriterCallEbos( *this, timer, state, wellState, cellData, substep ) );
+            }
+            else {
+                // just write the data to disk
+                writeTimeStepSerial( timer, state, wellState, cellData, substep );
+            }
+        }
+    }
+
+
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepSerial(const SimulatorTimerInterface& timer,
+                        const SimulationDataContainer& state,
+                        const WellState& wellState,
+                        const std::vector<data::CellData>& simProps,
+                        bool substep)
+    {
+        // ECL output
+        if ( eclWriter_ )
+        {
+            const auto& initConfig = eclipseState_->getInitConfig();
+            if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
+                std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
+            } else {
+                eclWriter_->writeTimeStep(timer.reportStepNum(),
+                                          substep,
+                                          timer.simulationTimeElapsed(),
+                                          simToSolution( state, phaseUsage_ ),
+                                          wellState.report(),
+                                          simProps);
+            }
+        }
+
+        // write backup file
+        if( backupfile_.is_open() )
+        {
+            int reportStep      = timer.reportStepNum();
+            int currentTimeStep = timer.currentStepNum();
+            if( (reportStep == currentTimeStep || // true for SimulatorTimer
+                 currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep
+                 timer.done() ) // true for AdaptiveSimulatorTimer at reportStep
+               && lastBackupReportStep_ != reportStep ) // only backup report step once
+            {
+                // store report step
+                lastBackupReportStep_ = reportStep;
+                // write resport step number
+                backupfile_.write( (const char *) &reportStep, sizeof(int) );
+
+                try {
+                    backupfile_ << state;
+
+                    const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
+                    backupfile_ << boWellState;
+                }
+                catch ( const std::bad_cast& e )
+                {
+                }
+
+                backupfile_ << std::flush;
+            }
+        } // end backup
+    }
+
+    void
+    BlackoilOutputWriterEbos::
+    restore(SimulatorTimerInterface& timer,
+            BlackoilState& state,
+            WellStateFullyImplicitBlackoil& wellState,
+            const std::string& filename,
+            const int desiredResportStep )
+    {
+        std::ifstream restorefile( filename.c_str() );
+        if( restorefile )
+        {
+            std::cout << "============================================================================"<<std::endl;
+            std::cout << "Restoring from ";
+            if( desiredResportStep < 0 ) {
+                std::cout << "last";
+            }
+            else {
+                std::cout << desiredResportStep;
+            }
+            std::cout << " report step! filename = " << filename << std::endl << std::endl;
+
+            int reportStep;
+            restorefile.read( (char *) &reportStep, sizeof(int) );
+
+            const int readReportStep = (desiredResportStep < 0) ?
+                std::numeric_limits<int>::max() : desiredResportStep;
+
+            while( reportStep <= readReportStep && ! timer.done() && restorefile )
+            {
+                restorefile >> state;
+                restorefile >> wellState;
+
+                // No per cell data is written for restore steps, but will be
+                // for subsequent steps, when we have started simulating
+                writeTimeStepWithoutCellProperties( timer, state, wellState );
+
+                // some output
+                std::cout << "Restored step " << timer.reportStepNum() << " at day "
+                          <<  unit::convert::to(timer.simulationTimeElapsed(),unit::day) << std::endl;
+
+                if( readReportStep == reportStep ) {
+                    break;
+                }
+
+                // if the stream is not valid anymore we just use the last state read
+                if( ! restorefile ) {
+                    std::cerr << "Reached EOF, using last state read!" << std::endl;
+                    break;
+                }
+
+                // try to read next report step
+                restorefile.read( (char *) &reportStep, sizeof(int) );
+
+                // if read failed, exit loop
+                if( ! restorefile ) {
+                    break;
+                }
+
+                // next step
+                timer.advance();
+
+                if( timer.reportStepNum() != reportStep ) {
+                    break;
+                }
+            }
+        }
+        else
+        {
+            std::cerr << "Warning: Couldn't open restore file '" << filename << "'" << std::endl;
+        }
+    }
+
+
+    bool BlackoilOutputWriterEbos::isRestart() const {
+        const auto& initconfig = eclipseState_->getInitConfig();
+        return initconfig.restartRequested();
+    }
+}

opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp

@@ -0,0 +1,500 @@
+/*
+  Copyright (c) 2014 SINTEF ICT, Applied Mathematics.
+  Copyright (c) 2015 IRIS AS
+
+  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_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED
+#include <opm/core/grid.h>
+#include <opm/core/simulator/SimulatorTimerInterface.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/Compat.hpp>
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/output/eclipse/EclipseReader.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/wells/DynamicListEconLimited.hpp>
+
+#include <opm/output/Cells.hpp>
+#include <opm/output/eclipse/EclipseWriter.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/ParallelDebugOutput.hpp>
+
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/ThreadHandle.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+
+
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <thread>
+
+#include <boost/filesystem.hpp>
+
+#ifdef HAVE_OPM_GRID
+#include <dune/grid/CpGrid.hpp>
+#endif
+namespace Opm
+{
+    class BlackoilState;
+
+    /** \brief Wrapper class for VTK, Matlab, and ECL output. */
+    class BlackoilOutputWriterEbos
+    {
+
+    public:
+        // constructor creating different sub writers
+        template <class Grid>
+        BlackoilOutputWriterEbos(const Grid& grid,
+                                 const parameter::ParameterGroup& param,
+                                 Opm::EclipseStateConstPtr eclipseState,
+                                 const Opm::PhaseUsage &phaseUsage,
+                                 const double* permeability );
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will extract the
+         *        requested output cell properties specified by the RPTRST keyword
+         *        and write these to file.
+         */
+        template<class Model>
+        void writeTimeStep(const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           const Model& physicalModel,
+                           bool substep = false);
+
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will write all
+         *        CellData in simProps to the file as well.
+         */
+        void writeTimeStepWithCellProperties(
+                           const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           const std::vector<data::CellData>& simProps,
+                           bool substep = false);
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will not write
+         *        any cell properties (e.g., those requested by RPTRST keyword)
+         */
+        void writeTimeStepWithoutCellProperties(
+                           const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           bool substep = false);
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection withS
+         *        visualization tools like ResInsight. This is the function which does
+         *        the actual write to file.
+         */
+        void writeTimeStepSerial(const SimulatorTimerInterface& timer,
+                                 const SimulationDataContainer& reservoirState,
+                                 const Opm::WellState& wellState,
+                                 const std::vector<data::CellData>& simProps,
+                                 bool substep);
+
+        /** \brief return output directory */
+        const std::string& outputDirectory() const { return outputDir_; }
+
+        /** \brief return true if output is enabled */
+        bool output () const { return output_; }
+
+        void restore(SimulatorTimerInterface& timer,
+                     BlackoilState& state,
+                     WellStateFullyImplicitBlackoil& wellState,
+                     const std::string& filename,
+                     const int desiredReportStep);
+
+
+        template <class Grid>
+        void initFromRestartFile(const PhaseUsage& phaseusage,
+                                 const double* permeability,
+                                 const Grid& grid,
+                                 SimulationDataContainer& simulatorstate,
+                                 WellStateFullyImplicitBlackoil& wellstate);
+
+        bool isRestart() const;
+
+    protected:
+        const bool output_;
+        std::unique_ptr< ParallelDebugOutputInterface > parallelOutput_;
+
+        // Parameters for output.
+        const std::string outputDir_;
+        const int output_interval_;
+
+        int lastBackupReportStep_;
+
+        std::ofstream backupfile_;
+        Opm::PhaseUsage phaseUsage_;
+        std::unique_ptr< EclipseWriter > eclWriter_;
+        EclipseStateConstPtr eclipseState_;
+
+        std::unique_ptr< ThreadHandle > asyncOutput_;
+    };
+
+
+    //////////////////////////////////////////////////////////////
+    //
+    //  Implementation
+    //
+    //////////////////////////////////////////////////////////////
+    template <class Grid>
+    inline
+    BlackoilOutputWriterEbos::
+    BlackoilOutputWriterEbos(const Grid& grid,
+                         const parameter::ParameterGroup& param,
+                         Opm::EclipseStateConstPtr eclipseState,
+                         const Opm::PhaseUsage &phaseUsage,
+                         const double* permeability )
+      : output_( param.getDefault("output", true) ),
+        parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability ) : 0 ),
+        outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ),
+        output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ),
+        lastBackupReportStep_( -1 ),
+        phaseUsage_( phaseUsage ),
+        eclWriter_( output_ && parallelOutput_->isIORank() &&
+                    param.getDefault("output_ecl", true) ?
+                    new EclipseWriter(eclipseState,UgGridHelpers::createEclipseGrid( grid , *eclipseState->getInputGrid()))
+                   : 0 ),
+        eclipseState_(eclipseState),
+        asyncOutput_()
+    {
+        // For output.
+        if (output_ && parallelOutput_->isIORank() ) {
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(outputDir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+
+            // create output thread if enabled and rank is I/O rank
+            // async output is enabled by default if pthread are enabled
+#if HAVE_PTHREAD
+            const bool asyncOutputDefault = false;
+#else
+            const bool asyncOutputDefault = false;
+#endif
+            if( param.getDefault("async_output", asyncOutputDefault ) )
+            {
+#if HAVE_PTHREAD
+                asyncOutput_.reset( new ThreadHandle() );
+#else
+                OPM_THROW(std::runtime_error,"Pthreads were not found, cannot enable async_output");
+#endif
+            }
+
+            std::string backupfilename = param.getDefault("backupfile", std::string("") );
+            if( ! backupfilename.empty() )
+            {
+                backupfile_.open( backupfilename.c_str() );
+            }
+        }
+    }
+
+
+    template <class Grid>
+    inline void
+    BlackoilOutputWriterEbos::
+    initFromRestartFile( const PhaseUsage& phaseusage,
+                         const double* permeability,
+                         const Grid& grid,
+                         SimulationDataContainer& simulatorstate,
+                         WellStateFullyImplicitBlackoil& wellstate)
+    {
+        // gives a dummy dynamic_list_econ_limited
+        DynamicListEconLimited dummy_list_econ_limited;
+        WellsManager wellsmanager(eclipseState_,
+                                  eclipseState_->getInitConfig().getRestartStep(),
+                                  Opm::UgGridHelpers::numCells(grid),
+                                  Opm::UgGridHelpers::globalCell(grid),
+                                  Opm::UgGridHelpers::cartDims(grid),
+                                  Opm::UgGridHelpers::dimensions(grid),
+                                  Opm::UgGridHelpers::cell2Faces(grid),
+                                  Opm::UgGridHelpers::beginFaceCentroids(grid),
+                                  permeability,
+                                  dummy_list_econ_limited);
+
+        const Wells* wells = wellsmanager.c_wells();
+        wellstate.resize(wells, simulatorstate); //Resize for restart step
+        auto restarted = Opm::init_from_restart_file(
+                                *eclipseState_,
+                                Opm::UgGridHelpers::numCells(grid) );
+
+        solutionToSim( restarted.first, phaseusage, simulatorstate );
+        wellsToState( restarted.second, wellstate );
+    }
+
+
+
+
+
+    namespace detail {
+        template<class Model>
+        std::vector<data::CellData> getCellDataEbos(
+                const Opm::PhaseUsage& phaseUsage,
+                const Model& model,
+                const RestartConfig& restartConfig,
+                const int reportStepNum)
+        {
+            typedef typename Model::FluidSystem FluidSystem;
+
+            std::vector<data::CellData> simProps;
+
+            //Get the value of each of the keys
+            std::map<std::string, int> outKeywords = restartConfig.getRestartKeywords(reportStepNum);
+            for (auto& keyValue : outKeywords) {
+                keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
+            }
+
+            //Get shorthands for water, oil, gas
+            const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
+            const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
+            const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
+
+            const auto& ebosModel = model.ebosSimulator().model();
+
+            // extract everything which can possibly be written to disk
+            int numCells = ebosModel.numGridDof();
+            std::vector<double> bWater(numCells);
+            std::vector<double> bOil(numCells);
+            std::vector<double> bGas(numCells);
+
+            std::vector<double> rhoWater(numCells);
+            std::vector<double> rhoOil(numCells);
+            std::vector<double> rhoGas(numCells);
+
+            std::vector<double> muWater(numCells);
+            std::vector<double> muOil(numCells);
+            std::vector<double> muGas(numCells);
+
+            std::vector<double> krWater(numCells);
+            std::vector<double> krOil(numCells);
+            std::vector<double> krGas(numCells);
+
+            std::vector<double> Rs(numCells);
+            std::vector<double> Rv(numCells);
+
+            for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+                const auto& intQuants = *ebosModel.cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
+                const auto& fs = intQuants.fluidState();
+
+                bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value;
+                bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value;
+                bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value;
+
+                rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value;
+                rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value;
+                rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value;
+
+                muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value;
+                muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value;
+                muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value;
+
+                krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value;
+                krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value;
+                krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value;
+
+                Rs[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
+                                                                      FluidSystem::oilPhaseIdx,
+                                                                      intQuants.pvtRegionIndex(),
+                                                                      /*maxOilSaturation=*/1.0).value;
+                Rv[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
+                                                                      FluidSystem::gasPhaseIdx,
+                                                                      intQuants.pvtRegionIndex(),
+                                                                      /*maxOilSaturation=*/1.0).value;
+            }
+
+            /**
+             * Formation volume factors for water, oil, gas
+             */
+            if (aqua_active && outKeywords["BW"] > 0) {
+                outKeywords["BW"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "1OVERBW",
+                        Opm::UnitSystem::measure::water_inverse_formation_volume_factor,
+                        std::move(bWater)});
+            }
+            if (liquid_active && outKeywords["BO"]  > 0) {
+                outKeywords["BO"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "1OVERBO",
+                        Opm::UnitSystem::measure::oil_inverse_formation_volume_factor,
+                        std::move(bOil)});
+            }
+            if (vapour_active && outKeywords["BG"] > 0) {
+                outKeywords["BG"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "1OVERBG",
+                        Opm::UnitSystem::measure::gas_inverse_formation_volume_factor,
+                        std::move(bGas)});
+            }
+
+            /**
+             * Densities for water, oil gas
+             */
+            if (outKeywords["DEN"] > 0) {
+                outKeywords["DEN"] = 0;
+                if (aqua_active) {
+                    simProps.emplace_back(data::CellData{
+                            "WAT_DEN",
+                            Opm::UnitSystem::measure::density,
+                            std::move(rhoWater)});
+                }
+                if (liquid_active) {
+                    simProps.emplace_back(data::CellData{
+                            "OIL_DEN",
+                            Opm::UnitSystem::measure::density,
+                            std::move(rhoOil)});
+                }
+                if (vapour_active) {
+                    simProps.emplace_back(data::CellData{
+                            "GAS_DEN",
+                            Opm::UnitSystem::measure::density,
+                            std::move(rhoGas)});
+                }
+            }
+
+            /**
+             * Viscosities for water, oil gas
+             */
+            if (outKeywords["VISC"] > 0) {
+                outKeywords["VISC"] = 0;
+                if (aqua_active) {
+                    simProps.emplace_back(data::CellData{
+                            "WAT_VISC",
+                            Opm::UnitSystem::measure::viscosity,
+                            std::move(muWater)});
+                }
+                if (liquid_active) {
+                    simProps.emplace_back(data::CellData{
+                            "OIL_VISC",
+                            Opm::UnitSystem::measure::viscosity,
+                            std::move(muOil)});
+                }
+                if (vapour_active) {
+                    simProps.emplace_back(data::CellData{
+                            "GAS_VISC",
+                            Opm::UnitSystem::measure::viscosity,
+                            std::move(muGas)});
+                }
+            }
+
+            /**
+             * Relative permeabilities for water, oil, gas
+             */
+            if (aqua_active && outKeywords["KRW"] > 0) {
+                outKeywords["KRW"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "WATKR",
+                        Opm::UnitSystem::measure::permeability,
+                        std::move(krWater)});
+            }
+            if (liquid_active && outKeywords["KRO"] > 0) {
+                outKeywords["KRO"] = 0;
+                simProps.emplace_back(data::CellData{
+                         "OILKR",
+                         Opm::UnitSystem::measure::permeability,
+                         std::move(krOil)});
+            }
+            if (vapour_active && outKeywords["KRG"] > 0) {
+                outKeywords["KRG"] = 0;
+                simProps.emplace_back(data::CellData{
+                          "GASKR",
+                          Opm::UnitSystem::measure::permeability,
+                          std::move(krGas)});
+            }
+
+            /**
+             * Vaporized and dissolved gas/oil ratio
+             */
+            if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
+                outKeywords["RSSAT"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "RSSAT",
+                        Opm::UnitSystem::measure::gas_oil_ratio,
+                        std::move(Rs)});
+            }
+            if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) {
+                outKeywords["RVSAT"] = 0;
+                simProps.emplace_back(data::CellData{
+                        "RVSAT",
+                        Opm::UnitSystem::measure::oil_gas_ratio,
+                        std::move(Rv)});
+            }
+
+
+            /**
+             * Bubble point and dew point pressures
+             */
+            if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) {
+                outKeywords["PBPD"] = 0;
+                Opm::OpmLog::warning("Bubble/dew point pressure output unsupported",
+                        "Writing bubble points and dew points (PBPD) to file is unsupported, "
+                        "as the simulator does not use these internally.");
+            }
+
+            //Warn for any unhandled keyword
+            for (auto& keyValue : outKeywords) {
+                if (keyValue.second > 0) {
+                    std::string logstring = "Keyword '";
+                    logstring.append(keyValue.first);
+                    logstring.append("' is unhandled for output to file.");
+                    Opm::OpmLog::warning("Unhandled output keyword", logstring);
+                }
+            }
+
+            return simProps;
+        }
+    }
+
+
+
+
+    template<class Model>
+    inline void
+    BlackoilOutputWriterEbos::
+    writeTimeStep(const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  const Model& physicalModel,
+                  bool substep)
+    {
+        const RestartConfig& restartConfig = eclipseState_->getRestartConfig();
+        const int reportStepNum = timer.reportStepNum();
+        std::vector<data::CellData> cellData = detail::getCellDataEbos( phaseUsage_, physicalModel, restartConfig, reportStepNum );
+        writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep);
+    }
+}
+#endif

opm/autodiff/StandardWellsDense.hpp

@@ -47,7 +47,7 @@
 #include <opm/autodiff/WellDensitySegmented.hpp>
 #include <opm/autodiff/BlackoilDetails.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
-
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 
 #include<dune/common/fmatrix.hh>
 #include<dune/istl/bcrsmatrix.hh>
@@ -305,7 +305,7 @@ namespace Opm {
             }
 
             void addRhs(BVector& x, Mat& jac) const {
-                assert(x.size() == rhs.size());
+                assert(x.size() == rhs_.size());
                 x += rhs_;
                 // jac = A + duneA
                 jac = matAdd( jac, duneA_ );
@@ -955,7 +955,6 @@ namespace Opm {
                                 dx_new_eigen(idx) = dx_new[w][flowPhaseToEbosCompIdx(p)];
                             }
                         }
-                        assert(dx.size() == total_residual_v.size());
                         updateWellState(dx_new_eigen.array(), well_state);
                         updateWellControls(well_state);
                         setWellVariables(well_state);
@@ -1995,7 +1994,6 @@ namespace Opm {
             }
 
             EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const {
-                assert(fluid_.numPhases() == 3);
                 const int nw = wells().number_of_wells;
                 if (phaseIdx == Water) {
                    return wellVariables_[nw + wellIdx];

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -225,7 +225,6 @@ namespace Opm
                                 }
                             }
                         }
-
                         // perfPressures
                         if( num_perf_old_well == num_perf_this_well )
                         {

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp

@@ -205,8 +205,8 @@ namespace {
          PolymerBlackoilState& x ,
          WellStateFullyImplicitBlackoilPolymer& xw)
     {
-        const double dt = timer.currentStepLength();
         const std::vector<double>& polymer_inflow = xw.polymerInflow();
+        const double dt = timer.currentStepLength();
 
         // Initial max concentration of this time step from PolymerBlackoilState.
         cmax_ = Eigen::Map<V>(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_));