Merge pull request #1289 from totto82/removeWellStateFullyImplicitBlackoilDense.hpp

Remove WellStateFullyImplicitBlackoilDense

CMakeLists_files.cmake

@@ -233,7 +233,6 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/TransportSolverTwophaseAd.hpp
   opm/autodiff/WellDensitySegmented.hpp
   opm/autodiff/WellStateFullyImplicitBlackoil.hpp
-  opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
   opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
   opm/autodiff/VFPProperties.hpp

opm/autodiff/BlackoilModelEbos.hpp

@@ -110,7 +110,7 @@ namespace Opm {
     public:
         // ---------  Types and enums  ---------
         typedef BlackoilState ReservoirState;
-        typedef WellStateFullyImplicitBlackoilDense WellState;
+        typedef WellStateFullyImplicitBlackoil WellState;
         typedef BlackoilModelParameters ModelParameters;
 
         typedef typename GET_PROP_TYPE(TypeTag, Simulator)         Simulator;

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -1098,7 +1098,8 @@ namespace Opm {
         if (has_solvent_ && is_miscible_ && b_eff_[0].size() == 0) {
             // A hack to avoid trouble for initial fluid in place, due to usage of b_eff_.
             WellState xw, xwdummy;
-            xw.init(&wells(), x, xwdummy);
+            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+            xw.init(&wells(), x, xwdummy, pu);
             SolutionState solstate = variableState(x, xw);
             computeEffectiveProperties(solstate);
         }

opm/autodiff/Compat.cpp

@@ -30,7 +30,6 @@
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
 #include <opm/output/data/Cells.hpp>
 #include <opm/output/data/Solution.hpp>
 #include <opm/output/data/Wells.hpp>
@@ -283,18 +282,4 @@ void wellsToState( const data::Wells& wells,
     }
 }
 
-
-
-
-
-
-void wellsToState( const data::Wells& wells,
-                   PhaseUsage phases,
-                   WellStateFullyImplicitBlackoilDense& state )
-{
-    // Set base class variables.
-    wellsToState(wells, phases, static_cast<WellStateFullyImplicitBlackoil&>(state));
-}
-
-
 } // namespace Opm

opm/autodiff/ParallelDebugOutput.hpp

@@ -621,7 +621,7 @@ namespace Opm
                                            std::unordered_set<std::string>());
 
                 const Wells* wells = wells_manager.c_wells();
-                globalWellState_.init(wells, *globalReservoirState_, globalWellState_ );
+                globalWellState_.init(wells, *globalReservoirState_, globalWellState_, phaseUsage_ );
                 globalCellData_->clear();
             }
 

opm/autodiff/SimulatorBase_impl.hpp

@@ -184,7 +184,7 @@ namespace Opm
                                        defunct_well_names_);
             const Wells* wells = wells_manager.c_wells();
             WellState well_state;
-            well_state.init(wells, state, prev_well_state);
+            well_state.init(wells, state, prev_well_state, props_.phaseUsage());
 
             // give the polymer and surfactant simulators the chance to do their stuff
             asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -26,7 +26,7 @@
 #include <opm/autodiff/NonlinearSolver.hpp>
 #include <opm/autodiff/BlackoilModelEbos.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/StandardWellsDense.hpp>
 #include <opm/autodiff/RateConverter.hpp>
 #include <opm/autodiff/SimFIBODetails.hpp>
@@ -60,7 +60,7 @@ public:
 
     typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
 
-    typedef WellStateFullyImplicitBlackoilDense WellState;
+    typedef WellStateFullyImplicitBlackoil WellState;
     typedef BlackoilState ReservoirState;
     typedef BlackoilOutputWriter OutputWriter;
     typedef BlackoilModelEbos<TypeTag> Model;

opm/autodiff/StandardWellsDense.hpp

@@ -45,7 +45,7 @@
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
 #include <opm/autodiff/BlackoilDetails.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/RateConverter.hpp>
 #include <opm/autodiff/WellInterface.hpp>
 #include <opm/autodiff/StandardWell.hpp>
@@ -65,7 +65,7 @@ namespace Opm {
         class StandardWellsDense {
         public:
             // ---------      Types      ---------
-            typedef WellStateFullyImplicitBlackoilDense WellState;
+            typedef WellStateFullyImplicitBlackoil WellState;
             typedef BlackoilModelParameters ModelParameters;
 
             typedef typename GET_PROP_TYPE(TypeTag, Grid)                Grid;

opm/autodiff/WellInterface.hpp

@@ -35,7 +35,7 @@
 #include <opm/autodiff/VFPInjProperties.hpp>
 #include <opm/autodiff/VFPProdProperties.hpp>
 #include <opm/autodiff/WellHelpers.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 
 #include <opm/simulators/WellSwitchingLogger.hpp>
@@ -61,7 +61,7 @@ namespace Opm
     {
     public:
 
-        using WellState = WellStateFullyImplicitBlackoilDense;
+        using WellState = WellStateFullyImplicitBlackoil;
 
         typedef BlackoilModelParameters ModelParameters;
         typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -54,16 +55,16 @@ namespace Opm
         using BaseType :: numPhases;
 
         template <class State, class PrevWellState>
-        void init(const Wells* wells, const State& state, const PrevWellState& prevState)
+        void init(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
         {
-            init(wells, state.pressure(), prevState);
+            init(wells, state.pressure(), prevState, pu);
         }
 
         /// 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 PrevWellState>
-        void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState)
+        void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
         {
             // call init on base class
             BaseType :: init(wells, cellPressures);
@@ -111,6 +112,9 @@ namespace Opm
 
             is_new_well_.resize(nw, true);
 
+            perfRateSolvent_.clear();
+            perfRateSolvent_.resize(nperf, 0.0);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if( ! prevState.wellMap().empty() )
@@ -141,18 +145,18 @@ namespace Opm
                         }
 
                         // perfPhaseRates
-                        int oldPerf_idx = (*it).second[ 1 ];
+                        const int oldPerf_idx_beg = (*it).second[ 1 ];
                         const int num_perf_old_well = (*it).second[ 2 ];
                         const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
                         // copy perforation rates when the number of perforations is equal,
                         // otherwise initialize perfphaserates to well rates divided by the number of perforations.
                         if( num_perf_old_well == num_perf_this_well )
                         {
-                            int oldPerf = oldPerf_idx *np;
+                            int old_perf_phase_idx = oldPerf_idx_beg *np;
-                            for (int perf = wells->well_connpos[ newIndex ]*np;
+                            for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
-                                 perf < wells->well_connpos[ newIndex + 1]*np; ++perf, ++oldPerf )
+                                 perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
                             {
-                                perfPhaseRates()[ perf ] = prevState.perfPhaseRates()[ oldPerf ];
+                                perfPhaseRates()[ perf_phase_idx ] = prevState.perfPhaseRates()[ old_perf_phase_idx ];
                             }
                         } else {
                             for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
@@ -161,15 +165,28 @@ namespace Opm
                                 }
                             }
                         }
-                        // perfPressures
+                        // perfPressures                        
                         if( num_perf_old_well == num_perf_this_well )
                         {
+                            int oldPerf_idx = oldPerf_idx_beg;
                             for (int perf = wells->well_connpos[ newIndex ];
                                  perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
                             {
                                 perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
                             }
                         }
+                        // perfSolventRates
+                        if (pu.has_solvent) {
+                            if( num_perf_old_well == num_perf_this_well )
+                            {
+                                int oldPerf_idx = oldPerf_idx_beg;
+                                for (int perf = wells->well_connpos[ newIndex ];
+                                     perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
+                                {
+                                    perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
+                                }
+                            }
+                        }
                     }
 
 
@@ -192,14 +209,9 @@ namespace Opm
         }
 
         template <class State>
-        void resize(const Wells* wells, const State& state ) {
+        void resize(const Wells* wells, const State& state, const PhaseUsage& pu) {
             const WellStateFullyImplicitBlackoil dummy_state{}; // Init with an empty previous state only resizes
-            init(wells, state, dummy_state) ;
+            init(wells, state, dummy_state, pu) ;
-        }
-
-        template <class State>
-        void resize(const Wells* wells, const State& state, const PhaseUsage& ) {
-            resize( wells, state );
         }
 
         /// One rate per phase and well connection.
@@ -232,6 +244,14 @@ namespace Opm
                 phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
             }
 
+            if (pu.has_solvent) {
+                // add solvent component
+                for( int w = 0; w < nw; ++w ) {
+                    using rt = data::Rates::opt;
+                    res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
+                }
+            }
+
             /* this is a reference or example on **how** to convert from
              * WellState to something understood by opm-output. it is intended
              * to be properly implemented and maintained as a part of
@@ -271,9 +291,24 @@ namespace Opm
             is_new_well_[w] = is_new_well;
         }
 
+
+        /// One rate pr well connection.
+        std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
+        const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
+
+        /// One rate pr well
+        double solventWellRate(const int w) const {
+            double solvent_well_rate = 0.0;
+            for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
+                solvent_well_rate += perfRateSolvent_[perf];
+            }
+            return solvent_well_rate;
+        }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
+        std::vector<double> perfRateSolvent_;
 
         // marking whether the well is just added
         // for newly added well, the current initialized rates from WellState

opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp

@@ -1,162 +0,0 @@
-/*
-  Copyright 2016 IRIS
-
-  This file is part of the Open Porous Media project (OPM).
-
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED
-#define OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED
-
-
-#include <opm/core/wells.h>
-#include <opm/core/well_controls.h>
-#include <opm/core/simulator/WellState.hpp>
-#include <opm/autodiff/BlackoilModelEnums.hpp>
-#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
-#include <opm/core/props/BlackoilPhases.hpp>
-#include <opm/common/ErrorMacros.hpp>
-#include <vector>
-#include <cassert>
-#include <string>
-#include <utility>
-#include <map>
-#include <algorithm>
-#include <array>
-#include <cmath>
-
-namespace Opm
-{
-
-    /// The state of a set of wells, tailored for use by the fully
-    /// implicit blackoil simulator.
-    class WellStateFullyImplicitBlackoilDense
-        : public WellStateFullyImplicitBlackoil
-    {
-        typedef WellStateFullyImplicitBlackoil  BaseType;
-    public:
-        typedef BaseType :: WellMapType WellMapType;
-
-        using BaseType :: wellRates;
-        using BaseType :: bhp;
-        using BaseType :: perfPress;
-        using BaseType :: wellMap;
-        using BaseType :: numWells;
-        using BaseType :: numPhases;
-        using BaseType :: perfPhaseRates;
-        using BaseType :: currentControls;
-
-        /// 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 PrevWellState>
-        void init(const Wells* wells, const std::vector<double>& cellPressures, const PrevWellState& prevState, const PhaseUsage& pu)
-        {
-
-            // call init on base class
-            BaseType :: init(wells, cellPressures, prevState);
-
-
-            const int nw = wells->number_of_wells;
-            if (nw == 0) {
-                return;
-            }
-            const int nperf = wells->well_connpos[nw];
-            perfRateSolvent_.clear();
-            perfRateSolvent_.resize(nperf, 0.0);
-
-            if (pu.has_solvent) {
-
-                // intialize wells that have been there before
-                // order may change so the mapping is based on the well name
-                if( ! prevState.wellMap().empty() )
-                {
-                    typedef typename WellMapType :: const_iterator const_iterator;
-                    const_iterator end = prevState.wellMap().end();
-                    for (int w = 0; w < nw; ++w) {
-                        std::string name( wells->name[ w ] );
-                        const_iterator it = prevState.wellMap().find( name );
-                        if( it != end )
-                        {
-                            const int newIndex = w;
-
-                            // perfSolventRates
-                            int oldPerf_idx = (*it).second[ 1 ];
-                            const int num_perf_old_well = (*it).second[ 2 ];
-                            const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
-                            if( num_perf_old_well == num_perf_this_well )
-                            {
-                                for (int perf = wells->well_connpos[ newIndex ];
-                                     perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
-                                {
-                                    perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-
-        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            
-            init(wells, state.pressure(), dummy_state, pu) ;
-        }
-
-        /// One rate pr well connection.
-        std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
-        const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
-
-        /// One rate pr well
-        double solventWellRate(const int w) const {
-            double solvent_well_rate = 0.0;
-            for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
-                solvent_well_rate += perfRateSolvent_[perf];
-            }
-            return solvent_well_rate;
-        }
-
-
-        data::Wells report(const PhaseUsage& pu) const override {
-            data::Wells res = BaseType::report(pu);
-            const int nw = WellState::numWells();
-            // If there are now wells numPhases throws a floating point
-            // exception.
-            if (nw == 0) {
-                return res;
-            }
-            if (pu.has_solvent) {
-                // add solvent component
-                for( int w = 0; w < nw; ++w ) {
-                    using rt = data::Rates::opt;
-                    res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
-                }
-            }
-
-            return res;
-        }
-
-
-    private:
-        std::vector<double> perfRateSolvent_;
-
-    };
-
-} // namespace Opm
-
-
-#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED