Merge pull request #1247 from totto82/removeState-PR

Remove reservoirState from BlackoilModelEbos

opm/autodiff/BlackoilModelEbos.hpp

@@ -115,23 +115,23 @@ namespace Opm {
         typedef typename GET_PROP_TYPE(TypeTag, SolutionVector)    SolutionVector ;
         typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables)  PrimaryVariables ;
         typedef typename GET_PROP_TYPE(TypeTag, FluidSystem)       FluidSystem;
-        typedef typename GET_PROP_TYPE(TypeTag, Indices)           BlackoilIndices;
+        typedef typename GET_PROP_TYPE(TypeTag, Indices)           Indices;
         typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)       MaterialLaw;
         typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
 
         typedef double Scalar;
-        static const int numEq = BlackoilIndices::numEq;
+        static const int numEq = Indices::numEq;
-        static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
+        static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
-        static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
+        static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
-        static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
+        static const int solventSaturationIdx = Indices::solventSaturationIdx;
-        static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
+        static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
 
         typedef Dune::FieldVector<Scalar, numEq >        VectorBlockType;
         typedef Dune::FieldMatrix<Scalar, numEq, numEq >        MatrixBlockType;
         typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
         typedef Dune::BlockVector<VectorBlockType>      BVector;
 
-        typedef ISTLSolver< MatrixBlockType, VectorBlockType, BlackoilIndices::pressureSwitchIdx >  ISTLSolverType;
+        typedef ISTLSolver< MatrixBlockType, VectorBlockType, Indices::pressureSwitchIdx >  ISTLSolverType;
         //typedef typename SolutionVector :: value_type            PrimaryVariables ;
 
         // For the conversion between the surface volume rate and resrevoir voidage rate
@@ -204,12 +204,16 @@ namespace Opm {
         /// \param[in, out] reservoir_state   reservoir state variables
         /// \param[in, out] well_state        well state variables
         void prepareStep(const SimulatorTimerInterface& /*timer*/,
-                         const ReservoirState& reservoir_state,
+                         const ReservoirState& /*reservoir_state*/,
                          const WellState& /* well_state */)
         {
             if ( wellModel().wellCollection()->havingVREPGroups() ) {
-                updateRateConverter(reservoir_state);
+                updateRateConverter();
             }
+
+            unsigned numDof = ebosSimulator_.model().numGridDof();
+            wasSwitched_.resize(numDof);
+            std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
         }
 
 
@@ -226,7 +230,7 @@ namespace Opm {
         SimulatorReport nonlinearIteration(const int iteration,
                                            const SimulatorTimerInterface& timer,
                                            NonlinearSolverType& nonlinear_solver,
-                                           ReservoirState& reservoir_state,
+                                           ReservoirState& /*reservoir_state*/,
                                            WellState& well_state)
         {
             SimulatorReport report;
@@ -245,7 +249,7 @@ namespace Opm {
             report.total_linearizations = 1;
 
             try {
-                report += assemble(timer, iteration, reservoir_state, well_state);
+                report += assemble(timer, iteration, well_state);
                 report.assemble_time += perfTimer.stop();
             }
             catch (...) {
@@ -298,30 +302,27 @@ namespace Opm {
                 perfTimer.reset();
                 perfTimer.start();
 
-                // Stabilize the nonlinear update.
+                if (param_.use_update_stabilization_) {
-                bool isOscillate = false;
+                    // Stabilize the nonlinear update.
-                bool isStagnate = false;
+                    bool isOscillate = false;
-                nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
+                    bool isStagnate = false;
-                if (isOscillate) {
+                    nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
-                    current_relaxation_ -= nonlinear_solver.relaxIncrement();
+                    if (isOscillate) {
-                    current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
+                        current_relaxation_ -= nonlinear_solver.relaxIncrement();
-                    if (terminalOutputEnabled()) {
+                        current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
-                        std::string msg = "    Oscillating behavior detected: Relaxation set to "
+                        if (terminalOutputEnabled()) {
-                            + std::to_string(current_relaxation_);
+                            std::string msg = "    Oscillating behavior detected: Relaxation set to "
-                        OpmLog::info(msg);
+                                    + std::to_string(current_relaxation_);
+                            OpmLog::info(msg);
+                        }
                     }
+                    nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
                 }
-                nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
 
                 // Apply the update, with considering model-dependent limitations and
                 // chopping of the update.
-                updateState(x,reservoir_state);
+                updateState(x,iteration);
                 wellModel().updateWellState(xw, well_state);
-                // if the solution is updated the solution needs to be comunicated to ebos
-                // and the cachedIntensiveQuantities needs to be updated.
-                convertInput( iteration, reservoir_state, ebosSimulator_ );
-                ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
-
                 report.update_time += perfTimer.stop();
             }
 
@@ -355,7 +356,6 @@ namespace Opm {
         /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
         SimulatorReport assemble(const SimulatorTimerInterface& timer,
                                  const int iterationIdx,
-                                 const ReservoirState& reservoir_state,
                                  WellState& well_state)
         {
             using namespace Opm::AutoDiffGrid;
@@ -364,11 +364,11 @@ namespace Opm {
 
             // when having VREP group control, update the rate converter based on reservoir state
             if ( wellModel().wellCollection()->havingVREPGroups() ) {
-                updateRateConverter(reservoir_state);
+                updateRateConverter();
             }
 
             // -------- Mass balance equations --------
-            assembleMassBalanceEq(timer, iterationIdx, reservoir_state);
+            assembleMassBalanceEq(timer, iterationIdx);
 
             // -------- Well equations ----------
             double dt = timer.currentStepLength();
@@ -389,44 +389,86 @@ namespace Opm {
             return report;
         }
 
-
+        // compute the "relative" change of the solution between time steps
-         /// \brief compute the relative change between to simulation states
+        template <class Dummy>
-        //  \return || u^n+1 - u^n || / || u^n+1 ||
+        double relativeChange(const Dummy&, const Dummy&) const
-        double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const
         {
-            std::vector< double > p0  ( previous.pressure() );
+            Scalar resultDelta = 0.0;
-            std::vector< double > sat0( previous.saturation() );
+            Scalar resultDenom = 0.0;
 
-            const std::size_t pSize = p0.size();
+            const auto& elemMapper = ebosSimulator_.model().elementMapper();
-            const std::size_t satSize = sat0.size();
+            const auto& gridView = ebosSimulator_.gridView();
+            auto elemIt = gridView.template begin</*codim=*/0>();
+            const auto& elemEndIt = gridView.template end</*codim=*/0>();
+            for (; elemIt != elemEndIt; ++elemIt) {
+                const auto& elem = *elemIt;
+                if (elem.partitionType() != Dune::InteriorEntity)
+                    continue;
 
-            // compute u^n - u^n+1
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
-            for( std::size_t i=0; i<pSize; ++i ) {
+		unsigned globalElemIdx = elemMapper.index(elem);
-                p0[ i ] -= current.pressure()[ i ];
+#else
+		unsigned globalElemIdx = elemMapper.map(elem);
+#endif
+                const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
+
+                Scalar pressureNew;
+		pressureNew = priVarsNew[Indices::pressureSwitchIdx];
+
+		Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
+                Scalar oilSaturationNew = 1.0;
+                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                    saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
+                    oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
+                }
+                
+                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
+                    saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
+                    oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];                     
+                }
+                                    
+                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                    saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
+                }
+
+                const auto& priVarsOld = ebosSimulator_.model().solution(/*timeIdx=*/1)[globalElemIdx];
+
+                Scalar pressureOld;
+                pressureOld = priVarsOld[Indices::pressureSwitchIdx];
+
+                Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
+                Scalar oilSaturationOld = 1.0;
+                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                    saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
+                    oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
+                }
+                
+                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
+                    saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
+                    oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
+                }
+                                    
+                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                    saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
+                }
+                
+                Scalar tmp = pressureNew - pressureOld;
+                resultDelta += tmp*tmp;
+                resultDenom += pressureNew*pressureNew;
+
+                for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
+                    Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
+		    resultDelta += tmp*tmp;
+                    resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
+                }
             }
 
-            for( std::size_t i=0; i<satSize; ++i ) {
+            resultDelta = gridView.comm().sum(resultDelta);
-                sat0[ i ] -= current.saturation()[ i ];
+            resultDenom = gridView.comm().sum(resultDenom);
-            }
 
-            // compute || u^n - u^n+1 ||
+	    if (resultDenom > 0.0)
-            const double stateOld  = detail::euclidianNormSquared( p0.begin(),   p0.end(), 1, istlSolver().parallelInformation() ) +
+	      return resultDelta/resultDenom;
-                detail::euclidianNormSquared( sat0.begin(), sat0.end(),
+	    return 0.0;
-                                              current.numPhases(),
-                                              istlSolver().parallelInformation() );
-
-            // compute || u^n+1 ||
-            const double stateNew  = detail::euclidianNormSquared( current.pressure().begin(),   current.pressure().end(), 1, istlSolver().parallelInformation() ) +
-                detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
-                                              current.numPhases(),
-                                              istlSolver().parallelInformation() );
-
-            if( stateNew > 0.0 ) {
-                return stateOld / stateNew ;
-            }
-            else {
-                return 0.0;
-            }
         }
 
 
@@ -567,25 +609,38 @@ namespace Opm {
         /// \param[in, out] reservoir_state   reservoir state variables
         /// \param[in, out] well_state        well state variables
         void updateState(const BVector& dx,
-                         ReservoirState& reservoir_state)
+                         const int iterationIdx)
         {
             using namespace Opm::AutoDiffGrid;
-            const int np = phaseUsage_.num_phases;
 
+            const auto& ebosProblem = ebosSimulator_.problem();
+
+            unsigned numSwitched = 0;
             ElementContext elemCtx( ebosSimulator_ );
             const auto& gridView = ebosSimulator_.gridView();
             const auto& elemEndIt = gridView.template end</*codim=*/0>();
+            SolutionVector& solution = ebosSimulator_.model().solution( 0 /* timeIdx */ );
+
+            // Store the initial solution.
+            if( iterationIdx == 0 )
+            {
+                ebosSimulator_.model().solution( 1 /* timeIdx */ ) = solution;
+            }
+
             for (auto elemIt = gridView.template begin</*codim=*/0>();
                  elemIt != elemEndIt;
                  ++elemIt)
             {
                 const auto& elem = *elemIt;
+
                 elemCtx.updatePrimaryStencil(elem);
                 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
                 const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+
+                PrimaryVariables& priVars = solution[ cell_idx ];
+
                 const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
-                //reservoir_state.pressure()[cell_idx] -= dp;
+                double& p = priVars[Indices::pressureSwitchIdx];
-                double& p = reservoir_state.pressure()[cell_idx];
                 const double& dp_rel_max = dpMaxRel();
                 const int sign_dp = dp > 0 ? 1: -1;
                 p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max);
@@ -601,169 +656,103 @@ namespace Opm {
                 double drs = 0.0;
                 double drv = 0.0;
 
-                double maxVal = 0.0;
+                // determine the saturation delta values
-                // water phase
+                if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
-                maxVal = std::max(std::abs(dsw),maxVal);
-                dso -= dsw;
-                // gas phase
-                switch (reservoir_state.hydroCarbonState()[cell_idx]) {
-                case HydroCarbonState::GasAndOil:
                     dsg = dxvar;
-                    break;
-                case HydroCarbonState::OilOnly:
-                    drs = dxvar;
-                    break;
-                case HydroCarbonState::GasOnly:
-                    dsg -= dsw;
-                    drv = dxvar;
-                    break;
-                default:
-                    OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << reservoir_state.hydroCarbonState()[cell_idx]);
                 }
-                dso -= dsg;
+                else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
+                    drs = dxvar;
+                }
+                else {
+                    assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
+                    drv = dxvar;
+                    dsg = 0.0;
+                }
 
                 // solvent
-                const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0;
+                const double dss = has_solvent_ ? dx[cell_idx][Indices::solventSaturationIdx] : 0.0;
-                dso -= dss;
 
                 // polymer
-                const double dc = has_polymer_ ? dx[cell_idx][BlackoilIndices::polymerConcentrationIdx] : 0.0;
+                const double dc = has_polymer_ ? dx[cell_idx][Indices::polymerConcentrationIdx] : 0.0;
 
-                // Appleyard chop process.
+                // oil
+                dso = - (dsw + dsg + dss);
+
+                // compute a scaling factor for the saturation update so that the maximum
+                // allowed change of saturations between iterations is not exceeded
+                double maxVal = 0.0;
+                maxVal = std::max(std::abs(dsw),maxVal);
                 maxVal = std::max(std::abs(dsg),maxVal);
-                maxVal = std::max(std::abs(dss),maxVal);
+                maxVal = std::max(std::abs(dso),maxVal);
+                maxVal = std::max(std::abs(dss),maxVal);                               
 
-                double step = dsMax()/maxVal;
+                double satScaleFactor = 1.0;
-                step = std::min(step, 1.0);
+                if (maxVal > dsMax()) {
-
+                    satScaleFactor = dsMax()/maxVal;
-                const Opm::PhaseUsage& pu = phaseUsage_;
-                if (active_[Water]) {
-                    double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
-                    sw -= step * dsw;
                 }
+
+                if (active_[Water]) {
+                    double& sw = priVars[Indices::waterSaturationIdx];
+                    sw -= satScaleFactor * dsw;
+                }
+
                 if (active_[Gas]) {
-                    double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
+                     if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
-                    sg -= step * dsg;
+                           double& sg = priVars[Indices::compositionSwitchIdx];
+                           sg -= satScaleFactor * dsg;
+                     }
                 }
 
                 if (has_solvent_) {
-                    double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
+                    double& ss = priVars[Indices::solventSaturationIdx];
-                    ss -= step * dss;
+                    ss -= satScaleFactor * dss;
                 }
-
                 if (has_polymer_) {
-                    double& c = reservoir_state.getCellData( reservoir_state.POLYMER )[cell_idx];
+                    double& c = priVars[Indices::polymerConcentrationIdx];
-                    c -= step * dc;
+                    c -= satScaleFactor * dc;
                     c = std::max(c, 0.0);
                 }
 
-                double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
+                // Update rs and rv
-                so -= step * dso;
-
-                // phase for when oil and gas
                 if (active_[Gas] && active_[Oil] ) {
-                    // const double drmaxrel = drMaxRel();
+                    unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
-                    // Update rs and rv
+                    const double drmaxrel = drMaxRel();
                     if (has_disgas_) {
-                        double& rs = reservoir_state.gasoilratio()[cell_idx];
+                        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
-                        rs -= drs;
+                            Scalar RsSat =
-                        rs = std::max(rs, 0.0);
+                                FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, 300.0, p);
+
+                            double& rs = priVars[Indices::compositionSwitchIdx];
+                            rs -= ((drs<0)?-1:1)*std::min(std::abs(drs), RsSat*drmaxrel);
+                            rs = std::max(rs, 0.0);
+                        }
 
                     }
                     if (has_vapoil_) {
-                        double& rv = reservoir_state.rv()[cell_idx];
+                        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
-                        rv -= drv;
+                            Scalar RvSat =
-                        rv = std::max(rv, 0.0);
+                                FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, 300.0, p);
-                    }
 
-                    // Sg is used as primal variable for water only cells.
+                            double& rv = priVars[Indices::compositionSwitchIdx];
-                    const double epsilon = 1e-4; //std::sqrt(std::numeric_limits<double>::epsilon());
+                            rv -= ((drv<0)?-1:1)*std::min(std::abs(drv), RvSat*drmaxrel);
-                    double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
+                            rv = std::max(rv, 0.0);
-                    double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
-                    double& rs = reservoir_state.gasoilratio()[cell_idx];
-                    double& rv = reservoir_state.rv()[cell_idx];
-
-                    // phase translation sg <-> rs
-                    const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
-                    const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
-                    const auto& fs = intQuants.fluidState();
-                    switch (hydroCarbonState) {
-                    case HydroCarbonState::GasAndOil: {
-
-                        // for the Gas and Oil case rs=rsSat and rv=rvSat
-                        rs = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
-                        // use gas pressure?
-                        rv = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
-
-                        if (sw > (1.0 - epsilon)) // water only i.e. do nothing
-                            break;
-
-                        if (sg <= 0.0 && has_disgas_) {
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
-                            sg = 0;
-                            so = 1.0 - sw;
-                            if (has_solvent_) {
-                                double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
-                                so -= ss;
-                            }
-                            rs *= (1-epsilon);
-                        } else if (so <= 0.0 && has_vapoil_) {
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv
-                            so = 0;
-                            sg = 1.0 - sw;
-                            if (has_solvent_) {
-                                double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
-                                sg -= ss;
-                            }
-                            rv *= (1-epsilon);
                         }
-                        break;
-                    }
-                    case HydroCarbonState::OilOnly: {
-                        if (sw > (1.0 - epsilon)) {
-                            // water only change to Sg
-                            rs = 0;
-                            rv = 0;
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
-                            //std::cout << "watonly rv -> sg" << cell_idx << std::endl;
-                            break;
-                        }
-
-                        const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
-                        if (rs > ( rsSat * (1+epsilon) ) ) {
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
-                            sg = epsilon;
-                            so -= epsilon;
-                            rs = rsSat;
-                        }
-                        break;
-                    }
-                    case HydroCarbonState::GasOnly: {
-                        if (sw > (1.0 - epsilon)) {
-                            // water only change to Sg
-                            rs = 0;
-                            rv = 0;
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
-                            //std::cout << "watonly rv -> sg" << cell_idx << std::endl;
-                            break;
-                        }
-
-                        const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
-                        if (rv > rvSat * (1+epsilon) ) {
-                            reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
-                            so = epsilon;
-                            rv = rvSat;
-                            sg -= epsilon;
-                        }
-                        break;
-                    }
-
-                    default:
-                        OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState);
                     }
                 }
+
+               // Add an epsilon to make it harder to switch back immediately after the primary variable was changed.
+                if (wasSwitched_[cell_idx])
+                    wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx, 1e-5);
+                else
+                    wasSwitched_[cell_idx] = priVars.adaptPrimaryVariables(ebosProblem, cell_idx);
+
+                if (wasSwitched_[cell_idx])
+                    ++numSwitched;
             }
 
+            // if the solution is updated the intensive Quantities need to be recalculated
+            ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
+
         }
 
         /// Return true if output to cout is wanted.
@@ -1568,114 +1557,47 @@ namespace Opm {
         bool localWellsActive() const { return well_model_.localWellsActive(); }
 
 
-        void convertInput( const int iterationIdx,
-                           const ReservoirState& reservoirState,
-                           Simulator& simulator ) const
-        {
-            SolutionVector& solution = simulator.model().solution( 0 /* timeIdx */ );
-            const Opm::PhaseUsage pu = phaseUsage_;
 
-            const int numCells = reservoirState.numCells();
-            const int numPhases = phaseUsage_.num_phases;
-            const auto& oilPressure = reservoirState.pressure();
-            const auto& saturations = reservoirState.saturation();
-            const auto& rs          = reservoirState.gasoilratio();
-            const auto& rv          = reservoirState.rv();
-            for( int cellIdx = 0; cellIdx<numCells; ++cellIdx )
-            {
-                // set non-switching primary variables
-                PrimaryVariables& cellPv = solution[ cellIdx ];
-                // set water saturation
-                cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
-
-                if (has_solvent_) {
-                    cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
-                }
-
-                if (has_polymer_) {
-                    cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
-                }
-
-
-                // set switching variable and interpretation
-                if (active_[Gas] ) {
-                    if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
-                    {
-                        cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
-                        cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
-                        cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
-                    }
-                    else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
-                    {
-                        // this case (-> gas only with vaporized oil in the gas) is
-                        // relatively expensive as it requires to compute the capillary
-                        // pressure in order to get the gas phase pressure. (the reason why
-                        // ebos uses the gas pressure here is that it makes the common case
-                        // of the primary variable switching code fast because to determine
-                        // whether the oil phase appears one needs to compute the Rv value
-                        // for the saturated gas phase and if this is not available as a
-                        // primary variable, it needs to be computed.) luckily for here, the
-                        // gas-only case is not too common, so the performance impact of this
-                        // is limited.
-                        typedef Opm::SimpleModularFluidState<double,
-                                /*numPhases=*/3,
-                                /*numComponents=*/3,
-                                FluidSystem,
-                                /*storePressure=*/false,
-                                /*storeTemperature=*/false,
-                                /*storeComposition=*/false,
-                                /*storeFugacity=*/false,
-                                /*storeSaturation=*/true,
-                                /*storeDensity=*/false,
-                                /*storeViscosity=*/false,
-                                /*storeEnthalpy=*/false> SatOnlyFluidState;
-                        SatOnlyFluidState fluidState;
-                        fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
-                        fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
-                        fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
-
-                        double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
-                        const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
-                        MaterialLaw::capillaryPressures(pC, matParams, fluidState);
-                        double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
-
-                        cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
-                        cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
-                        cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
-                    }
-                    else
-                    {
-                        assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
-                        cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
-                        cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
-                        cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
-                    }
-                } else {
-                    // for oil-water case oil pressure should be used as primary variable
-                    cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
-                }
-            }
-
-            if( iterationIdx == 0 )
-            {
-                simulator.model().solution( 1 /* timeIdx */ ) = solution;
-            }
-        }
 
     public:
+        int ebosCompToFlowPhaseIdx( const int compIdx ) const
+        {
+            assert(compIdx < 3);
+            const int compToPhase[ 3 ] = { Oil, Water, Gas };
+            return compToPhase[ compIdx ];
+        }
+
+        int flowToEbosPvIdx( const int flowPv ) const
+        {
+            const int flowToEbos[] = {
+                Indices::pressureSwitchIdx,
+                Indices::waterSaturationIdx,
+                Indices::compositionSwitchIdx,
+                Indices::solventSaturationIdx
+            };
+
+            if (flowPv > 2 )
+                return flowPv;
+
+            return flowToEbos[ flowPv ];
+        }
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const
         {
-            const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
+            const int phaseToComp[] = {
+                FluidSystem::waterCompIdx,
+                FluidSystem::oilCompIdx,
+                FluidSystem::gasCompIdx
+            };
+
             if (phaseIdx > 2 )
                 return phaseIdx;
+
             return phaseToComp[ phaseIdx ];
         }
 
-
     private:
 
-
         int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
         {
             assert(phaseIdx < 3);
@@ -1684,30 +1606,9 @@ namespace Opm {
         }
 
 
-        void updateRateConverter(const ReservoirState& reservoir_state)
+        void updateRateConverter()
         {
-            const int nw = numWells();
+            rate_converter_.defineState<ElementContext>(ebosSimulator_);
-            int global_number_wells = nw;
-
-#if HAVE_MPI
-            if ( istlSolver_->parallelInformation().type() == typeid(ParallelISTLInformation) )
-            {
-                const auto& info =
-                    boost::any_cast<const ParallelISTLInformation&>(istlSolver_->parallelInformation());
-                global_number_wells = info.communicator().sum(global_number_wells);
-                if ( global_number_wells )
-                {
-                    rate_converter_.defineState(reservoir_state, boost::any_cast<const ParallelISTLInformation&>(istlSolver_->parallelInformation()));
-                }
-            }
-            else
-#endif
-            {
-                if ( global_number_wells )
-                {
-                    rate_converter_.defineState(reservoir_state);
-                }
-            }
         }
 
 
@@ -1724,8 +1625,7 @@ namespace Opm {
 
     private:
         void assembleMassBalanceEq(const SimulatorTimerInterface& timer,
-                                   const int iterationIdx,
+                                   const int iterationIdx)
-                                   const ReservoirState& reservoirState)
         {
             ebosSimulator_.startNextEpisode( timer.currentStepLength() );
             ebosSimulator_.setEpisodeIndex( timer.reportStepNum() );
@@ -1755,9 +1655,9 @@ namespace Opm {
                 ebosSimulator_.problem().beginTimeStep();
             }
             // if the last time step failed we need to update the solution varables in ebos
-            // and recalculate the IntesiveQuantities. Also pass the solution initially.
+            // and recalculate the Intesive Quantities.
-            if ( (timer.lastStepFailed() || timer.reportStepNum()==0) && iterationIdx == 0  ) {
+            if ( timer.lastStepFailed() && iterationIdx == 0  ) {
-                convertInput( iterationIdx, reservoirState, ebosSimulator_ );
+                ebosSimulator_.model().solution( 0 /* timeIdx */ ) = ebosSimulator_.model().solution( 1 /* timeIdx */ );
                 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
             }
 
@@ -1780,6 +1680,7 @@ namespace Opm {
 
     public:
         bool isBeginReportStep_;
+        std::vector<bool> wasSwitched_;
     };
 } // namespace Opm
 

opm/autodiff/Compat.cpp

@@ -174,14 +174,17 @@ void solutionToSim( const data::Solution& sol,
     }
 
     if( sol.has( "RS" ) ) {
+        state.registerCellData("GASOILRATIO", 1);
         state.getCellData( "GASOILRATIO" ) = sol.data( "RS" );
     }
 
     if( sol.has( "RV" ) ) {
+        state.registerCellData("RV", 1);
         state.getCellData( "RV" ) = sol.data( "RV" );
     }
 
     if ( sol.has( "SSOL" ) ) {
+        state.registerCellData("SSOL", 1);
         state.getCellData("SSOL") = sol.data("SSOL");
     }
 

opm/autodiff/ParallelDebugOutput.hpp

@@ -345,8 +345,6 @@ namespace Opm
 
         class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface
         {
-            const SimulationDataContainer& localState_;
-            SimulationDataContainer& globalState_;
             const data::Solution& localCellData_;
             data::Solution& globalCellData_;
             const WellStateFullyImplicitBlackoil& localWellState_;
@@ -355,8 +353,7 @@ namespace Opm
             const IndexMapStorageType& indexMaps_;
 
         public:
-            PackUnPackSimulationDataContainer( const SimulationDataContainer& localState,
+            PackUnPackSimulationDataContainer( std::size_t numGlobalCells,
-                                               SimulationDataContainer& globalState,
                                                const data::Solution& localCellData,
                                                data::Solution& globalCellData,
                                                const WellStateFullyImplicitBlackoil& localWellState,
@@ -364,30 +361,20 @@ namespace Opm
                                                const IndexMapType& localIndexMap,
                                                const IndexMapStorageType& indexMaps,
                                                const bool isIORank )
-            : localState_( localState ),
+            : localCellData_( localCellData ),
-              globalState_( globalState ),
-              localCellData_( localCellData ),
               globalCellData_( globalCellData ),
               localWellState_( localWellState ),
               globalWellState_( globalWellState ),
               localIndexMap_( localIndexMap ),
               indexMaps_( indexMaps )
             {
+
                 if( isIORank )
                 {
-                    // add missing data to global state
-                    for (const auto& pair : localState.cellData()) {
-                        const std::string& key = pair.first;
-                        if (!globalState_.hasCellData( key )) {
-                            globalState_.registerCellData( key , localState.numCellDataComponents( key ));
-                        }
-                    }
-
                     // add missing data to global cell data
                     for (const auto& pair : localCellData_) {
                         const std::string& key = pair.first;
-                        std::size_t container_size = globalState_.numCells() *
+                        std::size_t container_size = numGlobalCells;
-                            pair.second.data.size() / localState_.numCells();
                         auto ret = globalCellData_.insert(key, pair.second.dim,
                                                 std::vector<double>(container_size),
                                                 pair.second.target);
@@ -413,13 +400,9 @@ namespace Opm
                 // write all cell data registered in local state
                 for (const auto& pair : localCellData_) {
                     const auto& data = pair.second.data;
-                    const size_t stride = data.size()/localState_.numCells();
 
-                    for( size_t i=0; i<stride; ++i )
+                    // write all data from local data to buffer
-                    {
+                    write( buffer, localIndexMap_, data);
-                        // write all data from local state to buffer
-                        write( buffer, localIndexMap_, data, i, stride );
-                    }
                 }
 
                 // write all data from local well state to buffer
@@ -428,23 +411,16 @@ namespace Opm
 
             void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
             {
-                // write all cell data registered in local state
+                // we loop over the data  as
-                // we loop over the data of the local state as
                 // its order governs the order the data got received.
                 for (auto& pair : localCellData_) {
                     const std::string& key = pair.first;
-
                     auto& data = globalCellData_.data(key);
-                    const size_t stride = data.size() / globalState_.numCells();
 
-                    for( size_t i=0; i<stride; ++i )
+                    //write all data from local cell data to buffer
-                    {
+                    read( buffer, indexMap, data);
-                        //write all data from local state to buffer
-                        read( buffer, indexMap, data, i, stride );
-                    }
                 }
 
-
                 // read well data from buffer
                 readWells( buffer );
             }
@@ -615,7 +591,7 @@ namespace Opm
         };
 
         // gather solution to rank 0 for EclipseWriter
-        bool collectToIORank( const SimulationDataContainer& localReservoirState,
+        bool collectToIORank( const SimulationDataContainer& /*localReservoirState*/,
                               const WellStateFullyImplicitBlackoil& localWellState,
                               const data::Solution& localCellData,
                               const int wellStateStepNumber )
@@ -649,7 +625,7 @@ namespace Opm
                 globalCellData_->clear();
             }
 
-            PackUnPackSimulationDataContainer packUnpack( localReservoirState, *globalReservoirState_,
+            PackUnPackSimulationDataContainer packUnpack( numCells(),
                                                           localCellData, *globalCellData_,
                                                           localWellState, globalWellState_,
                                                           localIndexMap_, indexMaps_,
@@ -663,7 +639,7 @@ namespace Opm
 #endif
             if( isIORank() )
             {
-                // Update values in the globalReservoirState
+                // copy values from globalCellData to globalReservoirState
                 const std::map<std::string, std::vector<double> > no_extra_data;
                 solutionToSim(*globalCellData_, no_extra_data, phaseUsage_, *globalReservoirState_);
             }

opm/autodiff/RateConverter.hpp

@@ -450,7 +450,66 @@ namespace Opm {
                 }   
                 calcRmax();
             }
+
             /**
+             * Compute average hydrocarbon pressure and maximum
+             * dissolution and evaporation at average hydrocarbon
+             * pressure in all regions in field.
+             *
+             * Fluid properties are evaluated at average hydrocarbon
+             * pressure for purpose of conversion from surface rate to
+             * reservoir voidage rate.
+             *
+             * \param[in] state Dynamic reservoir state.
+             * \param[in] any  The information and communication utilities
+             *                 about/of the parallelization. in any parallel
+             *                 it wraps a ParallelISTLInformation. Parameter
+             *                 is optional.
+             */
+            template <typename ElementContext, class EbosSimulator>
+            void defineState(const EbosSimulator& simulator)
+            {
+
+                //const int numCells = cellPvtIdx_.size();
+                //const Region region = std::vector<int>(numCells, 0);
+                auto& ra = attr_.attributes(0);
+                auto& p  = ra.pressure;
+                auto& T  = ra.temperature;
+                std::size_t n = 0;
+
+                ElementContext elemCtx( simulator );
+                const auto& gridView = simulator.gridView();
+
+                const auto& elemEndIt = gridView.template end</*codim=*/0>();
+                for (auto elemIt = gridView.template begin</*codim=*/0>();
+                     elemIt != elemEndIt;
+                     ++elemIt)
+                {
+
+                    const auto& elem = *elemIt;
+                    if (elem.partitionType() != Dune::InteriorEntity)
+                        continue;
+
+                    elemCtx.updatePrimaryStencil(elem);
+                    elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+                    const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+                    const auto& fs = intQuants.fluidState();
+
+                    p += fs.pressure(FluidSystem::oilPhaseIdx).value();
+                    T += fs.temperature(FluidSystem::oilPhaseIdx).value();
+                    n += 1;
+                }
+                p = gridView.comm().sum(p);
+                T = gridView.comm().sum(T);
+                n = gridView.comm().sum(n);
+
+                p /= n;
+                T /= n;
+
+                calcRmax();
+            }
+
+	    /**
              * Region identifier.
              *
              * Integral type.

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -55,6 +55,8 @@ public:
     typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables)  PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector)    SolutionVector ;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
 
     typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
 
@@ -113,7 +115,6 @@ public:
           defunct_well_names_( defunct_well_names ),
           is_parallel_run_( false )
     {
-
 #if HAVE_MPI
         if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
         {
@@ -143,11 +144,20 @@ public:
         extractLegacyPoreVolume_();
         extractLegacyDepth_();
 
+        // communicate the initial solution to ebos
+        if (timer.initialStep()) {
+            convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
+            ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
+        }
+
         if (output_writer_.isRestart()) {
             // This is a restart, populate WellState and ReservoirState state objects from restart file
             output_writer_.initFromRestartFile(phaseUsage_, grid(), state, prev_well_state, extra);
             initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
             initHysteresisParams(state);
+            // communicate the restart solution to ebos
+            convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
+            ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
         }
 
         // Create timers and file for writing timing info.
@@ -195,6 +205,11 @@ public:
                                     prev_well_state,
                                     restorefilename,
                                     desiredRestoreStep );
+            initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
+            initHysteresisParams(state);
+            // communicate the restart solution to ebos
+            convertInput(0, state, ebosSimulator_);
+            ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
         }
 
         DynamicListEconLimited dynamic_list_econ_limited;
@@ -239,13 +254,40 @@ public:
                                        defunct_well_names_ );
             const Wells* wells = wells_manager.c_wells();
             WellState well_state;
-            well_state.init(wells, state, prev_well_state, phaseUsage_);
+
+            // The well state initialize bhp with the cell pressure in the top cell.
+            // We must therefore provide it with updated cell pressures
+            size_t nc = Opm::UgGridHelpers::numCells(grid());
+            std::vector<double> cellPressures(nc, 0.0);
+            const auto& gridView = ebosSimulator_.gridManager().gridView();
+            ElementContext elemCtx(ebosSimulator_);
+            const auto& elemEndIt = gridView.template end</*codim=*/0>();
+            for (auto elemIt = gridView.template begin</*codim=*/0>();
+                 elemIt != elemEndIt;
+                 ++elemIt)
+            {
+                const auto& elem = *elemIt;
+                if (elem.partitionType() != Dune::InteriorEntity) {
+                    continue;
+                }
+
+                elemCtx.updatePrimaryStencil(elem);
+                elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+
+                const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+                const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+                const auto& fs = intQuants.fluidState();
+
+                const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
+                cellPressures[cellIdx] = p;
+            }
+            well_state.init(wells, cellPressures, prev_well_state, phaseUsage_);
 
             // give the polymer and surfactant simulators the chance to do their stuff
             handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
 
             // Compute reservoir volumes for RESV controls.
-            computeRESV(timer.currentStepNum(), wells, state, well_state);
+            computeRESV(timer.currentStepNum(), wells, well_state);
 
             // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
@@ -261,11 +303,8 @@ public:
 
             // Compute orignal fluid in place if this has not been done yet
             if (originalFluidInPlace.empty()) {
-                solver->model().convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
-                ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
-
                 originalFluidInPlace = solver->computeFluidInPlace(fipnum);
-                originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace, state);
+                originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace);
                 FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
                 FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
 
@@ -355,12 +394,19 @@ public:
                 stepReport.reportParam(tstep_os);
             }
 
+            // We don't need the reservoir state anymore. It is just passed around to avoid
+            // code duplication. Pass empty state instead.
+            if (timer.initialStep()) {
+                ReservoirState stateTrivial(0,0,0);
+                state = stateTrivial;
+            }
+
             // Increment timer, remember well state.
             ++timer;
 
             // Compute current fluid in place.
             currentFluidInPlace = solver->computeFluidInPlace(fipnum);
-            currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace, state);
+            currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace);
 
             const std::string version = moduleVersionName();
 
@@ -449,7 +495,6 @@ protected:
 
     void computeRESV(const std::size_t step,
                      const Wells* wells,
-                     const BlackoilState& x,
                      WellState& xw)
     {
         typedef SimFIBODetails::WellMap WellMap;
@@ -473,7 +518,7 @@ protected:
                 // to calculate averages over regions that might cross process
                 // borders. This needs to be done by all processes and therefore
                 // outside of the next if statement.
-                rateConverter_.defineState(x, boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()));
+                rateConverter_.template defineState<ElementContext>(ebosSimulator_);
             }
         }
         else
@@ -481,7 +526,7 @@ protected:
         {
             if ( global_number_resv_wells )
             {
-                rateConverter_.defineState(x);
+                rateConverter_.template defineState<ElementContext>(ebosSimulator_);
             }
         }
 
@@ -655,7 +700,7 @@ protected:
     }
 
 
-    std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip, const ReservoirState& /* state */)
+    std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip)
     {
         std::vector<double> totals(7,0.0);
         for (int i = 0; i < 5; ++i) {
@@ -850,6 +895,106 @@ protected:
         }
     }
 
+    // Used to convert initial Reservoirstate to primary variables in the SolutionVector
+    void convertInput( const int iterationIdx,
+                       const ReservoirState& reservoirState,
+                       Simulator& simulator ) const
+    {
+        SolutionVector& solution = simulator.model().solution( 0 /* timeIdx */ );
+        const Opm::PhaseUsage pu = phaseUsage_;
+
+        const std::vector<bool> active = detail::activePhases(pu);
+        bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
+        bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
+
+        const int numCells = reservoirState.numCells();
+        const int numPhases = phaseUsage_.num_phases;
+        const auto& oilPressure = reservoirState.pressure();
+        const auto& saturations = reservoirState.saturation();
+        const auto& rs          = reservoirState.gasoilratio();
+        const auto& rv          = reservoirState.rv();
+        for( int cellIdx = 0; cellIdx<numCells; ++cellIdx )
+        {
+            // set non-switching primary variables
+            PrimaryVariables& cellPv = solution[ cellIdx ];
+            // set water saturation
+            cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
+
+            if (has_solvent) {
+                cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
+            }
+
+            if (has_polymer) {
+                cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
+            }
+
+
+            // set switching variable and interpretation
+            if ( active[Gas] ) {
+                if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
+                {
+                    cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
+                    cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
+                    cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
+                }
+                else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
+                {
+                    // this case (-> gas only with vaporized oil in the gas) is
+                    // relatively expensive as it requires to compute the capillary
+                    // pressure in order to get the gas phase pressure. (the reason why
+                    // ebos uses the gas pressure here is that it makes the common case
+                    // of the primary variable switching code fast because to determine
+                    // whether the oil phase appears one needs to compute the Rv value
+                    // for the saturated gas phase and if this is not available as a
+                    // primary variable, it needs to be computed.) luckily for here, the
+                    // gas-only case is not too common, so the performance impact of this
+                    // is limited.
+                    typedef Opm::SimpleModularFluidState<double,
+                            /*numPhases=*/3,
+                            /*numComponents=*/3,
+                            FluidSystem,
+                            /*storePressure=*/false,
+                            /*storeTemperature=*/false,
+                            /*storeComposition=*/false,
+                            /*storeFugacity=*/false,
+                            /*storeSaturation=*/true,
+                            /*storeDensity=*/false,
+                            /*storeViscosity=*/false,
+                            /*storeEnthalpy=*/false> SatOnlyFluidState;
+                    SatOnlyFluidState fluidState;
+                    fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
+                    fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
+                    fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
+
+                    double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
+                    const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
+                    MaterialLaw::capillaryPressures(pC, matParams, fluidState);
+                    double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
+
+                    cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
+                    cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
+                    cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
+                }
+                else
+                {
+                    assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
+                    cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
+                    cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
+                    cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
+                }
+            } else {
+                // for oil-water case oil pressure should be used as primary variable
+                cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
+            }
+        }
+
+        // store the solution at the beginning of the time step
+        if( iterationIdx == 0 )
+        {
+            simulator.model().solution( 1 /* timeIdx */ ) = solution;
+        }
+    }
+
     RateConverterType createRateConverter_() {
         extractLegacyCellPvtRegionIndex_();
         RateConverterType rate_converter(phaseUsage_,

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -53,14 +53,20 @@ namespace Opm
         using BaseType :: numWells;
         using BaseType :: numPhases;
 
+        template <class State, class PrevWellState>
+        void init(const Wells* wells, const State& state, const PrevWellState& prevState)
+        {
+            init(wells, state.pressure(), prevState);
+        }
+
         /// Allocate and initialize if wells is non-null.  Also tries
         /// to give useful initial values to the bhp(), wellRates()
         /// and perfPhaseRates() fields, depending on controls
-        template <class State, class PrevState>
+        template <class PrevWellState>
-        void init(const Wells* wells, const State& state, const PrevState& prevState)
+        void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState)
         {
             // call init on base class
-            BaseType :: init(wells, state);
+            BaseType :: init(wells, cellPressures);
 
             // if there are no well, do nothing in init
             if (wells == 0) {
@@ -90,7 +96,7 @@ namespace Opm
                         for (int p = 0; p < np; ++p) {
                             perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
                         }
-                        perfPress()[perf] = state.pressure()[wells->well_cells[perf]];
+                        perfPress()[perf] = cellPressures[wells->well_cells[perf]];
                     }
                 }
             }

opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp

@@ -61,11 +61,12 @@ namespace Opm
         /// Allocate and initialize if wells is non-null.  Also tries
         /// to give useful initial values to the bhp(), wellRates()
         /// and perfPhaseRates() fields, depending on controls
-        template <class State, class PrevState>
+        template <class PrevWellState>
-        void init(const Wells* wells, const State& state, const PrevState& prevState, const PhaseUsage& pu)
+        void init(const Wells* wells, const std::vector<double>& cellPressures, const PrevWellState& prevState, const PhaseUsage& pu)
         {
+
             // call init on base class
-            BaseType :: init(wells, state, prevState);
+            BaseType :: init(wells, cellPressures, prevState);
 
 
             const int nw = wells->number_of_wells;
@@ -207,8 +208,8 @@ namespace Opm
 
         template <class State>
         void resize(const Wells* wells, const State& state, const PhaseUsage& pu ) {
-            const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes
+            const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes            
-            init(wells, state, dummy_state, pu) ;
+            init(wells, state.pressure(), dummy_state, pu) ;
         }