flow: let the wells be managed by EclProblem

opm/autodiff/BlackoilModelEbos.hpp

@@ -72,7 +72,6 @@ BEGIN_PROPERTIES
 
 NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem, FlowNonLinearSolver, FlowIstlSolver, FlowModelParameters, FlowTimeSteppingParameters));
 SET_STRING_PROP(EclFlowProblem, OutputDir, "");
-SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
 // default in flow is to formulate the equations in surface volumes
 SET_BOOL_PROP(EclFlowProblem, BlackoilConserveSurfaceVolume, true);
@@ -87,6 +86,8 @@ SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
 SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true);
 SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false);
 
+SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel<TypeTag>);
+
 END_PROPERTIES
 
 namespace Opm {
@@ -209,8 +210,6 @@ namespace Opm {
             wasSwitched_.resize(numDof);
             std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
 
-            wellModel().beginTimeStep(timer.reportStepNum(), timer.simulationTimeElapsed());
-
             if (param_.update_equations_scaling_) {
                 std::cout << "equation scaling not suported yet" << std::endl;
                 //updateEquationsScaling();
@@ -302,7 +301,7 @@ namespace Opm {
                 // handling well state update before oscillation treatment is a decision based
                 // on observation to avoid some big performance degeneration under some circumstances.
                 // there is no theorectical explanation which way is better for sure.
-                wellModel().recoverWellSolutionAndUpdateWellState(x);
+                wellModel().postSolve(x);
 
                 if (param_.use_update_stabilization_) {
                     // Stabilize the nonlinear update.
@@ -343,9 +342,7 @@ namespace Opm {
         /// \param[in] timer                  simulation timer
         void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer)
         {
-            wellModel().timeStepSucceeded(timer.simulationTimeElapsed());
             ebosSimulator_.problem().endTimeStep();
-
         }
 
         /// Assemble the residual and Jacobian of the nonlinear system.
@@ -361,32 +358,6 @@ namespace Opm {
             ebosSimulator_.model().linearizer().linearize();
             ebosSimulator_.problem().endIteration();
 
-            // -------- Current time step length ----------
-            const double dt = timer.currentStepLength();
-
-            // -------- Well equations ----------
-
-            try
-            {
-                // assembles the well equations and applies the wells to
-                // the reservoir equations as a source term.
-                wellModel().assemble(iterationIdx, dt);
-            }
-            catch ( const Dune::FMatrixError& )
-            {
-                OPM_THROW(Opm::NumericalIssue,"Error encounted when solving well equations");
-            }
-
-            auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
-            if (param_.matrix_add_well_contributions_) {
-                wellModel().addWellContributions(ebosJac);
-            }
-            if ( param_.preconditioner_add_well_contributions_ &&
-                 ! param_.matrix_add_well_contributions_ ) {
-                matrix_for_preconditioner_ .reset(new Mat(ebosJac));
-                wellModel().addWellContributions(*matrix_for_preconditioner_);
-            }
-
             return wellModel().lastReport();
         }
 
@@ -504,9 +475,6 @@ namespace Opm {
         /// r is the residual.
         void solveJacobianSystem(BVector& x) const
         {
-            const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
-            auto& ebosResid = ebosSimulator_.model().linearizer().residual();
-
             // J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
             // with the reservoir and D is the wells itself.
             // The full system is reduced to a number of cells X number of cells system via Schur complement
@@ -516,8 +484,8 @@ namespace Opm {
             // r = [r, r_well], where r is the residual and r_well the well residual.
             // r -= B^T * D^-1 r_well
 
-            // apply well residual to the residual.
+            auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
-            wellModel().apply(ebosResid);
+            auto& ebosResid = ebosSimulator_.model().linearizer().residual();
 
             // set initial guess
             x = 0.0;
@@ -1001,9 +969,9 @@ namespace Opm {
         const BlackoilWellModel<TypeTag>&
         wellModel() const { return well_model_; }
 
-        void beginReportStep()
+        void beginReportStep(bool isRestart)
         {
-            ebosSimulator_.problem().beginEpisode();
+            ebosSimulator_.problem().beginEpisode(isRestart);
         }
 
         void endReportStep()

opm/autodiff/BlackoilWellModel.hpp

@@ -24,6 +24,7 @@
 #ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
 #define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
 
+#include <ebos/eclproblem.hh>
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
@@ -60,12 +61,18 @@
 
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
+BEGIN_PROPERTIES
+
+NEW_PROP_TAG(EnableTerminalOutput);
+
+END_PROPERTIES
 
 namespace Opm {
 
         /// Class for handling the blackoil well model.
         template<typename TypeTag>
-        class BlackoilWellModel {
+        class BlackoilWellModel : public Ewoms::BaseAuxiliaryModule<TypeTag>
+        {
         public:
             // ---------      Types      ---------
             typedef WellStateFullyImplicitBlackoil WellState;
@@ -77,6 +84,11 @@ namespace Opm {
             typedef typename GET_PROP_TYPE(TypeTag, Indices)             Indices;
             typedef typename GET_PROP_TYPE(TypeTag, Simulator)           Simulator;
             typedef typename GET_PROP_TYPE(TypeTag, Scalar)              Scalar;
+            typedef typename GET_PROP_TYPE(TypeTag, RateVector)          RateVector;
+            typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector)      GlobalEqVector;
+            typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix)      JacobianMatrix;
+
+            typedef typename Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
 
             static const int numEq = Indices::numEq;
             static const int solventSaturationIdx = Indices::solventSaturationIdx;
@@ -101,51 +113,94 @@ namespace Opm {
             using RateConverterType = RateConverter::
                 SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
 
-            BlackoilWellModel(Simulator& ebosSimulator,
+            BlackoilWellModel(Simulator& ebosSimulator);
-                              const ModelParameters& param,
-                              const bool terminal_output);
 
-            void initFromRestartFile(const RestartValue& restartValues)
+            void init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule);
+
+            /////////////
+            // <eWoms auxiliary module stuff>
+            /////////////
+            unsigned numDofs() const
+            // No extra dofs are inserted for wells. (we use a Schur complement.)
+            { return 0; }
+
+            void addNeighbors(std::vector<NeighborSet>& neighbors) const;
+
+            void applyInitial()
+            {}
+
+            void linearize(JacobianMatrix& mat , GlobalEqVector& res);
+
+            void postSolve(GlobalEqVector& deltaX)
             {
-                // gives a dummy dynamic_list_econ_limited
+                recoverWellSolutionAndUpdateWellState(deltaX);
-                DynamicListEconLimited dummyListEconLimited;
-                const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
-                WellsManager wellsmanager(eclState(),
-                                          schedule(),
-                                          // The restart step value is used to identify wells present at the given
-                                          // time step. Wells that are added at the same time step as RESTART is initiated
-                                          // will not be present in a restart file. Use the previous time step to retrieve
-                                          // wells that have information written to the restart file.
-                                          std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
-                                          Opm::UgGridHelpers::numCells(grid()),
-                                          Opm::UgGridHelpers::globalCell(grid()),
-                                          Opm::UgGridHelpers::cartDims(grid()),
-                                          Opm::UgGridHelpers::dimensions(grid()),
-                                          Opm::UgGridHelpers::cell2Faces(grid()),
-                                          Opm::UgGridHelpers::beginFaceCentroids(grid()),
-                                          dummyListEconLimited,
-                                          grid().comm().size() > 1,
-                                          defunctWellNames);
-
-                const Wells* wells = wellsmanager.c_wells();
-
-                const int nw = wells->number_of_wells;
-                if (nw > 0) {
-                    auto phaseUsage = phaseUsageFromDeck(eclState());
-                    size_t numCells = Opm::UgGridHelpers::numCells(grid());
-                    well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step
-                    wellsToState(restartValues.wells, phaseUsage, well_state_);
-                    previous_well_state_ = well_state_;
-                }
             }
 
-            // compute the well fluxes and assemble them in to the reservoir equations as source terms
+            /////////////
-            // and in the well equations.
+            // </ eWoms auxiliary module stuff>
-            void assemble(const int iterationIdx,
+            /////////////
-                          const double dt);
 
-            // substract Binv(D)rw from r;
+            template <class Restarter>
-            void apply( BVector& r) const;
+            void deserialize(Restarter& res)
+            {
+                // TODO (?)
+            }
+
+            /*!
+             * \brief This method writes the complete state of the well
+             *        to the harddisk.
+             */
+            template <class Restarter>
+            void serialize(Restarter& res)
+            {
+                // TODO (?)
+            }
+
+            void beginEpisode(const Opm::EclipseState& eclState,
+                              const Opm::Schedule& schedule,
+                              bool isRestart)
+            {
+                size_t episodeIdx = ebosSimulator_.episodeIndex();
+                // beginEpisode in eclProblem advances the episode index
+                // we don't want this when we are at the beginning of an
+                // restart.
+                if (isRestart)
+                    episodeIdx -= 1;
+
+                beginReportStep(episodeIdx);
+            }
+
+            void beginTimeStep();
+
+            void beginIteration()
+            {
+                assemble(ebosSimulator_.model().newtonMethod().numIterations(),
+                         ebosSimulator_.timeStepSize());
+            }
+
+            void endIteration()
+            { }
+
+            void endTimeStep()
+            {
+                timeStepSucceeded(ebosSimulator_.time());
+            }
+
+            void endEpisode()
+            {
+                endReportStep();
+            }
+
+            template <class Context>
+            void computeTotalRatesForDof(RateVector& rate,
+                                         const Context& context,
+                                         unsigned spaceIdx,
+                                         unsigned timeIdx) const;
+
+            void initFromRestartFile(const RestartValue& restartValues);
+
+            Opm::data::Wells wellData() const
+            { return well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid())); }
 
             // subtract B*inv(D)*C * x from A*x
             void apply(const BVector& x, BVector& Ax) const;
@@ -153,10 +208,6 @@ namespace Opm {
             // apply well model with scaling of alpha
             void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
 
-            // using the solution x to recover the solution xw for wells and applying
-            // xw to update Well State
-            void recoverWellSolutionAndUpdateWellState(const BVector& x);
-
             // Check if well equations is converged.
             bool getWellConvergence(const std::vector<Scalar>& B_avg) const;
 
@@ -172,31 +223,13 @@ namespace Opm {
             // return the internal well state
             const WellState& wellState() const;
 
-            // only use this for restart.
+            const SimulatorReport& lastReport() const;
-            void setRestartWellState(const WellState& well_state);
-
-            // called at the beginning of a time step
-            void beginTimeStep(const int timeStepIdx,const double simulationTime);
-            // called at the end of a time step
-            void timeStepSucceeded(const double& simulationTime);
 
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
 
-            // called at the end of a report step
-            void endReportStep();
-
-            const SimulatorReport& lastReport() const;
-
-
-            void addWellContributions(Mat& mat)
-            {
-                for ( const auto& well: well_container_ ) {
-                        well->addWellContributions(mat);
-                }
-            }
-
         protected:
+
             void extractLegacyPressure_(std::vector<double>& cellPressure) const
             {
                 size_t nc = number_of_cells_;
@@ -234,6 +267,9 @@ namespace Opm {
             // a vector of all the wells.
             std::vector<WellInterfacePtr > well_container_;
 
+            // map from logically cartesian cell indices to compressed ones
+            std::vector<int> cartesian_to_compressed_;
+
             // create the well container
             std::vector<WellInterfacePtr > createWellContainer(const int time_step);
 
@@ -276,6 +312,21 @@ namespace Opm {
             const Schedule& schedule() const
             { return ebosSimulator_.vanguard().schedule(); }
 
+            // compute the well fluxes and assemble them in to the reservoir equations as source terms
+            // and in the well equations.
+            void assemble(const int iterationIdx,
+                          const double dt);
+
+            // called at the end of a time step
+            void timeStepSucceeded(const double& simulationTime);
+
+            // called at the end of a report step
+            void endReportStep();
+
+            // using the solution x to recover the solution xw for wells and applying
+            // xw to update Well State
+            void recoverWellSolutionAndUpdateWellState(const BVector& x);
+
             void updateWellControls();
 
             void updateGroupControls();
@@ -283,7 +334,7 @@ namespace Opm {
             // setting the well_solutions_ based on well_state.
             void updatePrimaryVariables();
 
-            void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
+            void setupCartesianToCompressed_(const int* global_cell, int number_of_cells);
 
             void computeRepRadiusPerfLength(const Grid& grid);
 
@@ -322,8 +373,7 @@ namespace Opm {
 
             void resetWellControlFromState() const;
 
-            void assembleWellEq(const double dt,
+            void assembleWellEq(const double dt);
-                                bool only_wells);
 
             // some preparation work, mostly related to group control and RESV,
             // at the beginning of each time step (Not report step)

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -4,16 +4,26 @@
 namespace Opm {
     template<typename TypeTag>
     BlackoilWellModel<TypeTag>::
-    BlackoilWellModel(Simulator& ebosSimulator,
+    BlackoilWellModel(Simulator& ebosSimulator)
-                      const ModelParameters& param,
-                      const bool terminal_output)
         : ebosSimulator_(ebosSimulator)
-        , param_(param)
-        , terminal_output_(terminal_output)
         , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
         , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
     {
-        const auto& eclState = ebosSimulator_.vanguard().eclState();
+        terminal_output_ = false;
+        if (ebosSimulator.gridView().comm().rank() == 0)
+            terminal_output_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
+    }
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule)
+    {
+        gravity_ = ebosSimulator_.problem().gravity()[2];
+
+        extractLegacyCellPvtRegionIndex_();
+        extractLegacyDepth_();
+
         phase_usage_ = phaseUsageFromDeck(eclState);
 
         const auto& gridView = ebosSimulator_.gridView();
@@ -28,6 +38,75 @@ namespace Opm {
         extractLegacyCellPvtRegionIndex_();
         extractLegacyDepth_();
         initial_step_ = true;
+
+        const auto& grid = ebosSimulator_.vanguard().grid();
+        const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
+        setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
+                                    cartDims[0]*cartDims[1]*cartDims[2]);
+
+        // add the eWoms auxiliary module for the wells to the list
+        ebosSimulator_.model().addAuxiliaryModule(this);
+    }
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    addNeighbors(std::vector<NeighborSet>& neighbors) const
+    {
+        if (!param_.matrix_add_well_contributions_) {
+            return;
+        }
+
+        // Create cartesian to compressed mapping
+        int last_time_step = schedule().getTimeMap().size() - 1;
+        const auto& schedule_wells = schedule().getWells();
+        const auto& cartesianSize = Opm::UgGridHelpers::cartDims(grid());
+
+        // initialize the additional cell connections introduced by wells.
+        for (const auto well : schedule_wells)
+        {
+            std::vector<int> wellCells;
+            // All possible connections of the well
+            const auto& connectionSet = well->getConnections(last_time_step);
+            wellCells.reserve(connectionSet.size());
+
+            for ( size_t c=0; c < connectionSet.size(); c++ )
+            {
+                const auto& connection = connectionSet.get(c);
+                int i = connection.getI();
+                int j = connection.getJ();
+                int k = connection.getK();
+                int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
+                int compressed_idx = cartesian_to_compressed_.at(cart_grid_idx);
+
+                if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
+                    wellCells.push_back(compressed_idx);
+                }
+            }
+
+            for (int cellIdx : wellCells) {
+                neighbors[cellIdx].insert(wellCells.begin(),
+                                          wellCells.end());
+            }
+        }
+    }
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    linearize(JacobianMatrix& mat , GlobalEqVector& res)
+    {
+        if (!localWellsActive())
+            return;
+
+        for (const auto& well: well_container_) {
+            if (param_.matrix_add_well_contributions_)
+                well->addWellContributions(mat);
+
+            // applying the well residual to reservoir residuals
+            // r = r - duneC_^T * invDuneD_ * resWell_
+            well->apply(res);
+        }
     }
 
 
@@ -126,15 +205,17 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    beginTimeStep(const int timeStepIdx, const double simulationTime) {
+    beginTimeStep() {
         well_state_ = previous_well_state_;
 
+        const int reportStepIdx = ebosSimulator_.episodeIndex();
+        const double simulationTime = ebosSimulator_.time();
 
         // test wells
-        wellTesting(timeStepIdx, simulationTime);
+        wellTesting(reportStepIdx, simulationTime);
 
         // create the well container
-        well_container_ = createWellContainer(timeStepIdx);
+        well_container_ = createWellContainer(reportStepIdx);
 
         // do the initialization for all the wells
         // TODO: to see whether we can postpone of the intialization of the well containers to
@@ -205,12 +286,6 @@ namespace Opm {
         }
     }
 
-    // only use this for restart.
-    template<typename TypeTag>
-    void
-    BlackoilWellModel<TypeTag>::
-    setRestartWellState(const WellState& well_state) { previous_well_state_ = well_state; }
-
     // called at the end of a report step
     template<typename TypeTag>
     void
@@ -238,6 +313,60 @@ namespace Opm {
         previous_well_state_ = well_state_;
     }
 
+
+    template<typename TypeTag>
+    template <class Context>
+    void
+    BlackoilWellModel<TypeTag>::
+    computeTotalRatesForDof(RateVector& rate,
+                            const Context& context,
+                            unsigned spaceIdx,
+                            unsigned timeIdx) const
+    {
+        rate = 0;
+        int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        for (const auto& well : well_container_)
+            well->addCellRates(rate, elemIdx);
+    }
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    initFromRestartFile(const RestartValue& restartValues)
+    {
+        // gives a dummy dynamic_list_econ_limited
+        DynamicListEconLimited dummyListEconLimited;
+        const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
+        WellsManager wellsmanager(eclState(),
+                                  schedule(),
+                                  // The restart step value is used to identify wells present at the given
+                                  // time step. Wells that are added at the same time step as RESTART is initiated
+                                  // will not be present in a restart file. Use the previous time step to retrieve
+                                  // wells that have information written to the restart file.
+                                  std::max(eclState().getInitConfig().getRestartStep() - 1, 0),
+                                  Opm::UgGridHelpers::numCells(grid()),
+                                  Opm::UgGridHelpers::globalCell(grid()),
+                                  Opm::UgGridHelpers::cartDims(grid()),
+                                  Opm::UgGridHelpers::dimensions(grid()),
+                                  Opm::UgGridHelpers::cell2Faces(grid()),
+                                  Opm::UgGridHelpers::beginFaceCentroids(grid()),
+                                  dummyListEconLimited,
+                                  grid().comm().size() > 1,
+                                  defunctWellNames);
+
+        const Wells* wells = wellsmanager.c_wells();
+
+        const int nw = wells->number_of_wells;
+        if (nw > 0) {
+            auto phaseUsage = phaseUsageFromDeck(eclState());
+            size_t numCells = Opm::UgGridHelpers::numCells(grid());
+            well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step
+            wellsToState(restartValues.wells, phaseUsage, well_state_);
+            previous_well_state_ = well_state_;
+        }
+        initial_step_ = false;
+    }
+
     template<typename TypeTag>
     std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
     BlackoilWellModel<TypeTag>::
@@ -401,7 +530,7 @@ namespace Opm {
             // basically, this is a more updated state from the solveWellEq based on fixed
             // reservoir state, will tihs be a better place to inialize the explict information?
         }
-        assembleWellEq(dt, false);
+        assembleWellEq(dt);
 
         last_report_.converged = true;
     }
@@ -413,32 +542,13 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    assembleWellEq(const double dt,
+    assembleWellEq(const double dt)
-                   bool only_wells)
     {
         for (auto& well : well_container_) {
-            well->assembleWellEq(ebosSimulator_, dt, well_state_, only_wells);
+            well->assembleWellEq(ebosSimulator_, dt, well_state_);
         }
     }
 
-    // applying the well residual to reservoir residuals
-    // r = r - duneC_^T * invDuneD_ * resWell_
-    template<typename TypeTag>
-    void
-    BlackoilWellModel<TypeTag>::
-    apply( BVector& r) const
-    {
-        if ( ! localWellsActive() ) {
-            return;
-        }
-
-        for (auto& well : well_container_) {
-            well->apply(r);
-        }
-    }
-
-
-
 
 
     // Ax = A x - C D^-1 B x
@@ -584,7 +694,7 @@ namespace Opm {
         int it  = 0;
         bool converged;
         do {
-            assembleWellEq(dt, true);
+            assembleWellEq(dt);
 
             converged = getWellConvergence(B_avg);
 
@@ -1071,16 +1181,17 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const
+    setupCartesianToCompressed_(const int* global_cell, int number_of_cartesian_cells)
     {
+        cartesian_to_compressed_.resize(number_of_cartesian_cells, -1);
         if (global_cell) {
-            for (int i = 0; i < number_of_cells; ++i) {
+            for (unsigned i = 0; i < number_of_cells_; ++i) {
-                cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
+                cartesian_to_compressed_[global_cell[i]] = i;
             }
         }
         else {
-            for (int i = 0; i < number_of_cells; ++i) {
+            for (unsigned i = 0; i < number_of_cells_; ++i) {
-                cartesian_to_compressed.insert(std::make_pair(i, i));
+                cartesian_to_compressed_[i] = i;
             }
         }
 
@@ -1091,16 +1202,8 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     computeRepRadiusPerfLength(const Grid& grid)
     {
-        // TODO, the function does not work for parallel running
-        // to be fixed later.
-        const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
-
-        std::map<int,int> cartesian_to_compressed;
-        setupCompressedToCartesian(global_cell, number_of_cells_,
-                                    cartesian_to_compressed);
-
         for (const auto& well : well_container_) {
-            well->computeRepRadiusPerfLength(grid, cartesian_to_compressed);
+            well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_);
         }
     }
 

opm/autodiff/MultisegmentWell.hpp

@@ -33,6 +33,7 @@ namespace Opm
     {
     public:
         typedef WellInterface<TypeTag> Base;
+        typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
 
         using typename Base::WellState;
         using typename Base::Simulator;
@@ -107,44 +108,43 @@ namespace Opm
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
-                          const int num_cells);
+                          const int num_cells) override;
 
 
-        virtual void initPrimaryVariablesEvaluation() const;
+        virtual void initPrimaryVariablesEvaluation() const override;
 
-        virtual void assembleWellEq(Simulator& ebosSimulator,
+        virtual void assembleWellEq(const Simulator& ebosSimulator,
                                     const double dt,
-                                    WellState& well_state,
+                                    WellState& well_state) override;
-                                    bool only_wells);
 
         /// updating the well state based the control mode specified with current
         // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(WellState& well_state) const;
+        virtual void updateWellStateWithTarget(WellState& well_state) const override;
 
         /// check whether the well equations get converged for this well
-        virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
+        virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
 
         /// Ax = Ax - C D^-1 B x
-        virtual void apply(const BVector& x, BVector& Ax) const;
+        virtual void apply(const BVector& x, BVector& Ax) const override;
         /// r = r - C D^-1 Rw
-        virtual void apply(BVector& r) const;
+        virtual void apply(BVector& r) const override;
 
         /// using the solution x to recover the solution xw for wells and applying
         /// xw to update Well State
         virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
-                                                           WellState& well_state) const;
+                                                           WellState& well_state) const override;
 
         /// computing the well potentials for group control
         virtual void computeWellPotentials(const Simulator& ebosSimulator,
                                            const WellState& well_state,
-                                           std::vector<double>& well_potentials);
+                                           std::vector<double>& well_potentials) override;
 
-        virtual void updatePrimaryVariables(const WellState& well_state) const;
+        virtual void updatePrimaryVariables(const WellState& well_state) const override;
 
-        virtual void solveEqAndUpdateWellState(WellState& well_state); // const?
+        virtual void solveEqAndUpdateWellState(WellState& well_state) override; // const?
 
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
-                                                 const WellState& well_state); // should be const?
+                                                 const WellState& well_state) override; // should be const?
 
         /// number of segments for this well
         /// int number_of_segments_;
@@ -152,6 +152,17 @@ namespace Opm
 
         int numberOfPerforations() const;
 
+        void addCellRates(RateVector& rates, int cellIdx) const override
+        {
+            for (int perfIdx = 0; perfIdx < number_of_perforations_; ++perfIdx) {
+                if (Base::cells()[perfIdx] == cellIdx) {
+                    for (int i = 0; i < RateVector::dimension; ++i) {
+                        rates[i] += connectionRates_[perfIdx][i];
+                    }
+                }
+            }
+        }
+
     protected:
         int number_segments_;
 
@@ -253,6 +264,8 @@ namespace Opm
 
         std::vector<double> segment_depth_diffs_;
 
+        std::vector<RateVector> connectionRates_;
+
         void initMatrixAndVectors(const int num_cells) const;
 
         // protected functions
@@ -337,14 +350,13 @@ namespace Opm
         bool accelerationalPressureLossConsidered() const;
 
         // TODO: try to make ebosSimulator const, as it should be
-        void iterateWellEquations(Simulator& ebosSimulator,
+        void iterateWellEquations(const Simulator& ebosSimulator,
                                   const double dt,
                                   WellState& well_state);
 
-        void assembleWellEqWithoutIteration(Simulator& ebosSimulator,
+        void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                             const double dt,
-                                            WellState& well_state,
+                                            WellState& well_state);
-                                            bool only_wells);
     };
 
 }

opm/autodiff/MultisegmentWell_impl.hpp

@@ -52,6 +52,10 @@ namespace Opm
         if (has_polymer) {
             OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
         }
+
+	if (Base::has_energy) {
+            OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
+	}
         // since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
         // for other facilities needed but not available from parser, we need to process them here
 
@@ -110,6 +114,8 @@ namespace Opm
     {
         Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
 
+        connectionRates_.resize(number_of_perforations_);
+
         // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
         // for MultisegmentWell, it is much more complicated.
         // It can be specified directly, it can be calculated from the segment depth,
@@ -228,10 +234,9 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    assembleWellEq(Simulator& ebosSimulator,
+    assembleWellEq(const Simulator& ebosSimulator,
                    const double dt,
-                   WellState& well_state,
+                   WellState& well_state)
-                   bool only_wells)
     {
 
         const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
@@ -239,7 +244,7 @@ namespace Opm
             iterateWellEquations(ebosSimulator, dt, well_state);
         }
 
-        assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
+        assembleWellEqWithoutIteration(ebosSimulator, dt, well_state);
     }
 
 
@@ -1719,7 +1724,7 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    iterateWellEquations(Simulator& ebosSimulator,
+    iterateWellEquations(const Simulator& ebosSimulator,
                          const double dt,
                          WellState& well_state)
     {
@@ -1732,7 +1737,7 @@ namespace Opm
         int it = 0;
         for (; it < max_iter_number; ++it) {
 
-            assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true);
+            assembleWellEqWithoutIteration(ebosSimulator, dt, well_state);
 
             const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
 
@@ -1762,19 +1767,16 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    assembleWellEqWithoutIteration(Simulator& ebosSimulator,
+    assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                    const double dt,
-                                   WellState& well_state,
+                                   WellState& well_state)
-                                   bool only_wells)
     {
         // calculate the fluid properties needed.
         computeSegmentFluidProperties(ebosSimulator);
 
         // clear all entries
-        if (!only_wells) {
+        duneB_ = 0.0;
-            duneB_ = 0.0;
+        duneC_ = 0.0;
-            duneC_ = 0.0;
-        }
 
         duneD_ = 0.0;
         resWell_ = 0.0;
@@ -1784,9 +1786,6 @@ namespace Opm
         //
         // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
 
-        auto& ebosJac = ebosSimulator.model().linearizer().matrix();
-        auto& ebosResid = ebosSimulator.model().linearizer().residual();
-
         const bool allow_cf = getAllowCrossFlow();
 
         const int nseg = numberOfSegments();
@@ -1835,31 +1834,24 @@ namespace Opm
                     // the cq_s entering mass balance equations need to consider the efficiency factors.
                     const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
 
-                    if (!only_wells) {
+                    connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
-                        // subtract sum of component fluxes in the reservoir equation.
-                        // need to consider the efficiency factor
-                        ebosResid[cell_idx][comp_idx] -= cq_s_effective.value();
-                    }
 
                     // subtract sum of phase fluxes in the well equations.
                     resWell_[seg][comp_idx] -= cq_s_effective.value();
 
                     // assemble the jacobians
                     for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-                        if (!only_wells) {
+
-                            // also need to consider the efficiency factor when manipulating the jacobians.
+                        // also need to consider the efficiency factor when manipulating the jacobians.
-                            duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
+                        duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
-                        }
+                        
                         // the index name for the D should be eq_idx / pv_idx
                         duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
                     }
 
                     for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
-                        if (!only_wells) {
+                        // also need to consider the efficiency factor when manipulating the jacobians.
-                            // also need to consider the efficiency factor when manipulating the jacobians.
+                        duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);   
-                            ebosJac[cell_idx][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
-                            duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
-                        }
                     }
                 }
                 // TODO: we should save the perforation pressure and preforation rates?

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -188,18 +188,13 @@ public:
         SimulatorReport report;
         SimulatorReport stepReport;
 
-        WellModel wellModel(ebosSimulator_, modelParam_, terminalOutput_);
         if (isRestart()) {
-            wellModel.initFromRestartFile(*restartValues);
+            wellModel_().initFromRestartFile(*restartValues);
         }
 
-        if (modelParam_.matrix_add_well_contributions_ ||
+        // beginReportStep(...) wants to know when we are at the
-             modelParam_.preconditioner_add_well_contributions_)
+        // beginning of a restart
-        {
+        bool firstRestartStep = isRestart();
-            ebosSimulator_.model().clearAuxiliaryModules();
-            wellAuxMod_.reset(new WellConnectionAuxiliaryModule<TypeTag>(schedule(), grid()));
-            ebosSimulator_.model().addAuxiliaryModule(wellAuxMod_.get());
-        }
 
         // Main simulation loop.
         while (!timer.done()) {
@@ -210,30 +205,25 @@ public:
                 OpmLog::debug(ss.str());
             }
 
-            // Run a multiple steps of the solver depending on the time step control.
-            solverTimer.start();
-            //ebosSimulator_.setEpisodeIndex(timer.reportStepNum());
-            wellModel.beginReportStep(timer.currentStepNum());
-
-            auto solver = createSolver(wellModel);
-
             // write the inital state at the report stage
             if (timer.initialStep()) {
                 Dune::Timer perfTimer;
                 perfTimer.start();
 
-                // No per cell data is written for initial step, but will be
+                wellModel_().beginReportStep(timer.currentStepNum());
-                // for subsequent steps, when we have started simulating
+                ebosSimulator_.problem().writeOutput(false);
-                auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()));
-                ebosSimulator_.problem().writeOutput(localWellData,
-                                                     timer.simulationTimeElapsed(),
-                                                     /*isSubstep=*/false,
-                                                     totalTimer.secsSinceStart(),
-                                                     /*nextStepSize=*/-1.0);
 
                 report.output_write_time += perfTimer.stop();
             }
 
+            // Run a multiple steps of the solver depending on the time step control.
+            solverTimer.start();
+
+            auto solver = createSolver(wellModel_());
+
+            solver->model().beginReportStep(firstRestartStep);
+            firstRestartStep = false;
+
             if (terminalOutput_) {
                 std::ostringstream stepMsg;
                 boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
@@ -246,8 +236,6 @@ public:
                 OpmLog::info(stepMsg.str());
             }
 
-            solver->model().beginReportStep();
-
             // If sub stepping is enabled allow the solver to sub cycle
             // in case the report steps are too large for the solver to converge
             //
@@ -282,7 +270,6 @@ public:
             }
 
             solver->model().endReportStep();
-            wellModel.endReportStep();
 
             // take time that was used to solve system for this reportStep
             solverTimer.stop();
@@ -305,13 +292,8 @@ public:
             Dune::Timer perfTimer;
             perfTimer.start();
             const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
-
+            ebosSimulator_.problem().setNextTimeStepSize(nextstep);
-            auto localWellData = wellModel.wellState().report(phaseUsage_, Opm::UgGridHelpers::globalCell(grid()));
+            ebosSimulator_.problem().writeOutput(false);
-            ebosSimulator_.problem().writeOutput(localWellData,
-                                                 timer.simulationTimeElapsed(),
-                                                 /*isSubstep=*/false,
-                                                 totalTimer.secsSinceStart(),
-                                                 nextstep);
             report.output_write_time += perfTimer.stop();
 
             if (terminalOutput_) {
@@ -379,6 +361,12 @@ protected:
         return initconfig.restartRequested();
     }
 
+    WellModel& wellModel_()
+    { return ebosSimulator_.problem().wellModel(); }
+
+    const WellModel& wellModel_() const
+    { return ebosSimulator_.problem().wellModel(); }
+
     // Data.
     Simulator& ebosSimulator_;
     std::unique_ptr<WellConnectionAuxiliaryModule<TypeTag>> wellAuxMod_;

opm/autodiff/StandardWell.hpp

@@ -37,6 +37,8 @@ namespace Opm
 
     public:
         typedef WellInterface<TypeTag> Base;
+        typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+
         // TODO: some functions working with AD variables handles only with values (double) without
         // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value.
         // And also, it can also be beneficial to make these functions hanle different types of AD variables.
@@ -129,52 +131,63 @@ namespace Opm
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
                           const double gravity_arg,
-                          const int num_cells);
+                          const int num_cells) override;
 
 
-        virtual void initPrimaryVariablesEvaluation() const;
+        virtual void initPrimaryVariablesEvaluation() const override;
 
-        virtual void assembleWellEq(Simulator& ebosSimulator,
+        virtual void assembleWellEq(const Simulator& ebosSimulator,
                                     const double dt,
-                                    WellState& well_state,
+                                    WellState& well_state) override;
-                                    bool only_wells);
 
         /// updating the well state based the control mode specified with current
         // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(WellState& well_state) const;
+        virtual void updateWellStateWithTarget(WellState& well_state) const override;
 
         /// check whether the well equations get converged for this well
-        virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
+        virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
 
         /// Ax = Ax - C D^-1 B x
-        virtual void apply(const BVector& x, BVector& Ax) const;
+        virtual void apply(const BVector& x, BVector& Ax) const override;
         /// r = r - C D^-1 Rw
-        virtual void apply(BVector& r) const;
+        virtual void apply(BVector& r) const override;
 
         /// using the solution x to recover the solution xw for wells and applying
         /// xw to update Well State
         virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
-                                                           WellState& well_state) const;
+                                                           WellState& well_state) const override;
 
         /// computing the well potentials for group control
         virtual void computeWellPotentials(const Simulator& ebosSimulator,
                                            const WellState& well_state,
-                                           std::vector<double>& well_potentials) /* const */;
+                                           std::vector<double>& well_potentials) /* const */ override;
 
-        virtual void updatePrimaryVariables(const WellState& well_state) const;
+        virtual void updatePrimaryVariables(const WellState& well_state) const override;
 
-        virtual void solveEqAndUpdateWellState(WellState& well_state);
+        virtual void solveEqAndUpdateWellState(WellState& well_state) override;
 
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
-                                                 const WellState& well_state); // should be const?
+                                                 const WellState& well_state) override; // should be const?
 
-        virtual void  addWellContributions(Mat& mat) const;
+        virtual void  addWellContributions(Mat& mat) const override;
 
         /// \brief Wether the Jacobian will also have well contributions in it.
-        virtual bool jacobianContainsWellContributions() const
+        virtual bool jacobianContainsWellContributions() const override
         {
             return param_.matrix_add_well_contributions_;
         }
+
+        void addCellRates(RateVector& rates, int cellIdx) const override
+        {
+            for (int perfIdx = 0; perfIdx < number_of_perforations_; ++perfIdx) {
+                if (Base::cells()[perfIdx] == cellIdx) {
+                    for (int i = 0; i < RateVector::dimension; ++i) {
+                        rates[i] += connectionRates_[perfIdx][i];
+                    }
+                }
+            }
+        }
+
     protected:
 
         // protected functions from the Base class
@@ -228,6 +241,8 @@ namespace Opm
         // diagonal matrix for the well
         DiagMatWell invDuneD_;
 
+        std::vector<RateVector> connectionRates_;
+
         // several vector used in the matrix calculation
         mutable BVectorWell Bx_;
         mutable BVectorWell invDrw_;

opm/autodiff/StandardWell_impl.hpp

@@ -58,6 +58,8 @@ namespace Opm
     {
         Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);
 
+        connectionRates_.resize(number_of_perforations_);
+
         perf_depth_.resize(number_of_perforations_, 0.);
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
             const int cell_idx = well_cells_[perf];
@@ -241,6 +243,7 @@ namespace Opm
     StandardWell<TypeTag>::
     wellSurfaceVolumeFraction(const int compIdx) const
     {
+
         EvalWell sum_volume_fraction_scaled = 0.;
         for (int idx = 0; idx < num_components_; ++idx) {
             sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
@@ -432,24 +435,18 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    assembleWellEq(Simulator& ebosSimulator,
+    assembleWellEq(const Simulator& ebosSimulator,
                    const double dt,
-                   WellState& well_state,
+                   WellState& well_state)
-                   bool only_wells)
     {
         const int np = number_of_phases_;
 
         // clear all entries
-        if (!only_wells) {
+        duneB_ = 0.0;
-            duneB_ = 0.0;
+        duneC_ = 0.0;
-            duneC_ = 0.0;
-        }
         invDuneD_ = 0.0;
         resWell_ = 0.0;
 
-        auto& ebosJac = ebosSimulator.model().linearizer().matrix();
-        auto& ebosResid = ebosSimulator.model().linearizer().residual();
-
         // TODO: it probably can be static member for StandardWell
         const double volume = 0.002831684659200; // 0.1 cu ft;
 
@@ -483,34 +480,28 @@ namespace Opm
                 well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
             }
 
+            if (has_energy) {
+                connectionRates_[perf][contiEnergyEqIdx] = 0.0;
+            }
+
             for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
                 // the cq_s entering mass balance equations need to consider the efficiency factors.
                 const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_;
 
-                if (!only_wells) {
+                connectionRates_[perf][componentIdx] = Base::restrictEval(cq_s_effective);
-                    // subtract sum of component fluxes in the reservoir equation.
-                    // need to consider the efficiency factor
-                    ebosResid[cell_idx][componentIdx] -= cq_s_effective.value();
-                }
 
                 // subtract sum of phase fluxes in the well equations.
                 resWell_[0][componentIdx] -= cq_s_effective.value();
 
                 // assemble the jacobians
                 for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
-                    if (!only_wells) {
+                    // also need to consider the efficiency factor when manipulating the jacobians.
-                        // also need to consider the efficiency factor when manipulating the jacobians.
+                    duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
-                        duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
-                    }
                     invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq);
                 }
 
                 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
-                    if (!only_wells) {
+                    duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
-                        // also need to consider the efficiency factor when manipulating the jacobians.
-                        ebosJac[cell_idx][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
-                        duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
-                    }
                 }
 
                 // Store the perforation phase flux for later usage.
@@ -574,15 +565,7 @@ namespace Opm
                     }
                     // compute the thermal flux
                     cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx));
-		    // scale the flux by the scaling factor for the energy equation
+                    connectionRates_[perf][contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
-                    cq_r_thermal *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor);
-
-                    if (!only_wells) {
-                        for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
-                            ebosJac[cell_idx][cell_idx][contiEnergyEqIdx][pvIdx] -= cq_r_thermal.derivative(pvIdx);
-                        }
-                        ebosResid[cell_idx][contiEnergyEqIdx] -= cq_r_thermal.value();
-                    }
                 }
             }
 
@@ -595,12 +578,7 @@ namespace Opm
                 } else {
                     cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
                 }
-                if (!only_wells) {
+                connectionRates_[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
-                    for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
-                        ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
-                    }
-                    ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
-                }
             }
 
             // Store the perforation pressure for later usage.
@@ -631,9 +609,16 @@ namespace Opm
         assembleControlEq();
 
         // do the local inversion of D.
-        // we do this manually with invertMatrix to always get our
+        try
-        // specializations in for 3x3 and 4x4 matrices.
+        {
-        Dune::ISTLUtility::invertMatrix(invDuneD_[0][0]);
+            Dune::ISTLUtility::invertMatrix(invDuneD_[0][0]);
+        }
+        catch( ... )
+        {
+            OPM_THROW(Opm::NumericalIssue,"Error when inverting local well equations for well " + name());
+        }
+
+
     }
 
 

opm/autodiff/WellInterface.hpp

@@ -82,6 +82,7 @@ namespace Opm
         typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
         typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
         typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+        typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
 
         static const int numEq = Indices::numEq;
         typedef double Scalar;
@@ -121,7 +122,7 @@ namespace Opm
         const int indexOfWell() const;
 
         /// Well cells.
-        const std::vector<int>& cells() {return well_cells_; }
+        const std::vector<int>& cells() const {return well_cells_; }
 
         /// Well type, INJECTOR or PRODUCER.
         WellType wellType() const;
@@ -142,20 +143,19 @@ namespace Opm
 
         virtual void solveEqAndUpdateWellState(WellState& well_state) = 0;
 
-        virtual void assembleWellEq(Simulator& ebosSimulator,
+        virtual void assembleWellEq(const Simulator& ebosSimulator,
                                     const double dt,
-                                    WellState& well_state,
+                                    WellState& well_state) = 0;
-                                    bool only_wells) = 0;
 
         void updateWellTestState(const WellState& well_state,
                                  const double& simulationTime,
                                  const bool& writeMessageToOPMLog,
-                                 WellTestState& wellTestState) const;
+                                 WellTestState& wellTestState
-
+                                 ) const;
 
         void setWellEfficiencyFactor(const double efficiency_factor);
 
-        void computeRepRadiusPerfLength(const Grid& grid, const std::map<int, int>& cartesian_to_compressed);
+        void computeRepRadiusPerfLength(const Grid& grid, const std::vector<int>& cartesian_to_compressed);
 
         /// using the solution x to recover the solution xw for wells and applying
         /// xw to update Well State
@@ -196,6 +196,20 @@ namespace Opm
         virtual void addWellContributions(Mat&) const
         {}
 
+        virtual void addCellRates(RateVector& rates, int cellIdx) const
+        {}
+
+        template <class EvalWell>
+        Eval restrictEval(const EvalWell& in) const
+        {
+            Eval out = 0.0;
+            out.setValue(in.value());
+            for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
+                out.setDerivative(eqIdx, in.derivative(eqIdx));
+            }
+            return out;
+        }
+
         void closeCompletions(WellTestState& wellTestState);
 
         const Well* wellEcl() const;

opm/autodiff/WellInterface_impl.hpp

@@ -875,7 +875,7 @@ namespace Opm
     void
     WellInterface<TypeTag>::
     computeRepRadiusPerfLength(const Grid& grid,
-                               const std::map<int, int>& cartesian_to_compressed)
+                               const std::vector<int>& cartesian_to_compressed)
     {
         const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
         auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
@@ -902,12 +902,12 @@ namespace Opm
 
                 const int* cpgdim = cart_dims;
                 const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
-                const std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
+                const int cell = cartesian_to_compressed[cart_grid_indx];
-                if (cgit == cartesian_to_compressed.end()) {
+
+                if (cell < 0) {
                     OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
                               << k << " not found in grid (well = " << name() << ')');
                 }
-                const int cell = cgit->second;
 
                 {
                     double radius = connection.rw();
@@ -1023,7 +1023,7 @@ namespace Opm
         bool converged;
         WellState well_state0 = well_state;
         do {
-            assembleWellEq(ebosSimulator, dt, well_state, true);
+            assembleWellEq(ebosSimulator, dt, well_state);
 
             auto report = getWellConvergence(B_avg);
             converged = report.converged();
@@ -1047,6 +1047,7 @@ namespace Opm
             if ( terminal_output ) {
                 OpmLog::debug("WellTest: Well equation for well" +name() + " solution failed in getting converged with " + std::to_string(it) + " iterations");
             }
+            well_state = well_state0;
         }
     }
 

opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp

@@ -203,7 +203,6 @@ namespace Opm {
 
             auto& ebosSimulator = solver.model().ebosSimulator();
             auto& ebosProblem = ebosSimulator.problem();
-            auto phaseUsage = Opm::phaseUsageFromDeck(ebosSimulator.vanguard().eclState());
 
             // create adaptive step timer with previously used sub step size
             AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
@@ -316,13 +315,7 @@ namespace Opm {
                         Opm::time::StopWatch perfTimer;
                         perfTimer.start();
 
-                        // The writeOutput expects a local data::solution vector and a local data::well vector.
+                        ebosProblem.writeOutput(/*isSubStep=*/true);
-                        auto localWellData = solver.model().wellModel().wellState().report(phaseUsage, Opm::UgGridHelpers::globalCell(ebosSimulator.vanguard().grid()));
-                        ebosProblem.writeOutput(localWellData,
-                                                substepTimer.simulationTimeElapsed(),
-                                                /*isSubstep=*/true,
-                                                substepReport.total_time,
-                                                /*nextStepSize=*/-1.0);
 
                         report.output_write_time += perfTimer.secsSinceStart();
                     }
@@ -361,6 +354,7 @@ namespace Opm {
 
                     ++restarts;
                 }
+                ebosProblem.setNextTimeStepSize(substepTimer.currentStepLength());
             }