Add polymer option to flow_ebos

No extra equation is added for polymer in the well equation.

Seperate executables are added for polymer: flow_ebos_polymer
and solvent: flow_ebos_solvent

Tested and verified on the test cases in polymer_test_suite

This PR should not effect the performance and results of the blackoil
simulator

CMakeLists_files.cmake

@@ -112,6 +112,8 @@ list (APPEND EXAMPLE_SOURCE_FILES
   examples/flow_legacy.cpp
   examples/flow_sequential.cpp
   examples/flow_ebos.cpp
+  examples/flow_ebos_solvent.cpp
+  examples/flow_ebos_polymer.cpp
   examples/flow_multisegment.cpp
   examples/flow_solvent.cpp
   examples/sim_2p_incomp.cpp
@@ -132,6 +134,8 @@ list (APPEND PROGRAM_SOURCE_FILES
   examples/sim_2p_incomp_ad.cpp
   examples/sim_2p_comp_reorder.cpp
   examples/flow_ebos.cpp
+  examples/flow_ebos_solvent.cpp
+  examples/flow_ebos_polymer.cpp
   examples/flow_legacy.cpp
   examples/flow_sequential.cpp
   examples/flow_solvent.cpp

examples/flow_ebos.cpp

@@ -34,6 +34,6 @@
 // ----------------- Main program -----------------
 int main(int argc, char** argv)
 {
-    Opm::FlowMainEbos mainfunc;
+    Opm::FlowMainEbos<TTAG(EclFlowProblem)> mainfunc;
     return mainfunc.execute(argc, argv);
 }

examples/flow_ebos_polymer.cpp

@@ -0,0 +1,42 @@
+/*
+  Copyright 2017 IRIS AS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
+#include <dune/grid/CpGrid.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/autodiff/FlowMainEbos.hpp>
+
+namespace Ewoms {
+namespace Properties {
+NEW_TYPE_TAG(EclFlowPolymerProblem, INHERITS_FROM(EclFlowProblem));
+SET_BOOL_PROP(EclFlowPolymerProblem, EnablePolymer, true);
+}}
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    Opm::FlowMainEbos<TTAG(EclFlowPolymerProblem)> mainfunc;
+    return mainfunc.execute(argc, argv);
+}

examples/flow_ebos_solvent.cpp

@@ -0,0 +1,42 @@
+/*
+  Copyright 2017 IRIS AS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
+#include <dune/grid/CpGrid.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/autodiff/FlowMainEbos.hpp>
+
+namespace Ewoms {
+namespace Properties {
+NEW_TYPE_TAG(EclFlowSolventProblem, INHERITS_FROM(EclFlowProblem));
+SET_BOOL_PROP(EclFlowSolventProblem, EnableSolvent, true);
+}}
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    Opm::FlowMainEbos<TTAG(EclFlowSolventProblem)> mainfunc;
+    return mainfunc.execute(argc, argv);
+}

opm/autodiff/BlackoilModelEbos.hpp

@@ -83,7 +83,6 @@ NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
 SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
 SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
 SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true);
-SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
 
 // SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy
 // code for fluid and satfunc handling gets fully retired.
@@ -97,6 +96,7 @@ namespace Opm {
     /// where gas can be dissolved in oil and vice versa. It
     /// uses an industry-standard TPFA discretization with per-phase
     /// upwind weighting of mobilities.
+    template <class TypeTag>
     class BlackoilModelEbos
     {
     public:
@@ -105,7 +105,6 @@ namespace Opm {
         typedef WellStateFullyImplicitBlackoilDense WellState;
         typedef BlackoilModelParameters ModelParameters;
 
-        typedef typename TTAG(EclFlowProblem) TypeTag;
         typedef typename GET_PROP_TYPE(TypeTag, Simulator)         Simulator;
         typedef typename GET_PROP_TYPE(TypeTag, Grid)              Grid;
         typedef typename GET_PROP_TYPE(TypeTag, ElementContext)    ElementContext;
@@ -118,7 +117,11 @@ namespace Opm {
 
         typedef double Scalar;
         static const int numEq = BlackoilIndices::numEq;
-        static const int solventCompIdx = 3; //TODO get this from ebos
+        static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
+        static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
+        static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
+        static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
+
         typedef Dune::FieldVector<Scalar, numEq >        VectorBlockType;
         typedef Dune::FieldMatrix<Scalar, numEq, numEq >        MatrixBlockType;
         typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
@@ -149,7 +152,8 @@ namespace Opm {
                           const ModelParameters& param,
                           const StandardWellsDense<TypeTag>& well_model,
                           const NewtonIterationBlackoilInterface& linsolver,
-                          const bool terminal_output)
+                          const bool terminal_output
+                          )
         : ebosSimulator_(ebosSimulator)
         , grid_(ebosSimulator_.gridManager().grid())
         , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
@@ -161,8 +165,9 @@ namespace Opm {
         , has_disgas_(FluidSystem::enableDissolvedGas())
         , has_vapoil_(FluidSystem::enableVaporizedOil())
         , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
+        , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
         , param_( param )
-        , well_model_ (well_model)
+        , well_model_ (well_model)        
         , terminal_output_ (terminal_output)
         , rate_converter_(phaseUsage_, ebosSimulator_.problem().pvtRegionArray().empty()?nullptr:ebosSimulator_.problem().pvtRegionArray().data(), AutoDiffGrid::numCells(grid_), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
         , current_relaxation_(1.0)
@@ -615,6 +620,9 @@ namespace Opm {
                 const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0;
                 dso -= dss;
 
+                // polymer
+                const double dc = has_polymer_ ? dx[cell_idx][BlackoilIndices::polymerConcentrationIdx] : 0.0;
+
                 // Appleyard chop process.
                 maxVal = std::max(std::abs(dsg),maxVal);
                 maxVal = std::max(std::abs(dss),maxVal);
@@ -637,6 +645,12 @@ namespace Opm {
                     ss -= step * dss;
                 }
 
+                if (has_polymer_) {
+                    double& c = reservoir_state.getCellData( reservoir_state.POLYMER )[cell_idx];
+                    c -= step * dc;
+                    c = std::max(c, 0.0);
+                }
+
                 double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
                 so -= step * dso;
 
@@ -867,10 +881,16 @@ namespace Opm {
                 }
 
                 if ( has_solvent_ ) {
-                    B_avg[ solventCompIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
+                    B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
-                    const auto R2 = ebosResid[cell_idx][solventCompIdx];
+                    const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
-                    R_sum[ solventCompIdx ] += R2;
+                    R_sum[ contiSolventEqIdx ] += R2;
-                    maxCoeff[ solventCompIdx ] = std::max( maxCoeff[ solventCompIdx ], std::abs( R2 ) / pvValue );
+                    maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
+                }
+                if (has_polymer_ ) {
+                     B_avg[ contiPolymerEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
+                    const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
+                    R_sum[ contiPolymerEqIdx ] += R2;
+                    maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
                 }
 
             }
@@ -938,7 +958,11 @@ namespace Opm {
                         key[ phaseIdx ] = std::toupper( phaseName.front() );
                     }
                     if (has_solvent_) {
-                        key[ solventCompIdx ] = "S";
+                        key[ solventSaturationIdx ] = "S";
+                    }
+
+                    if (has_polymer_) {
+                        key[ polymerConcentrationIdx ] = "P";
                     }
 
                     for (int compIdx = 0; compIdx < numComp; ++compIdx) {
@@ -1004,6 +1028,9 @@ namespace Opm {
             if (has_solvent_)
                 numComp ++;
 
+            if (has_polymer_)
+                numComp ++;
+
             return numComp;
         }
 
@@ -1450,6 +1477,7 @@ namespace Opm {
         const bool has_disgas_;
         const bool has_vapoil_;
         const bool has_solvent_;
+        const bool has_polymer_;
 
         ModelParameters                 param_;
         SimulatorReport failureReport_;
@@ -1513,6 +1541,11 @@ namespace Opm {
                     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_ )
@@ -1588,18 +1621,23 @@ namespace Opm {
 
         int flowToEbosPvIdx( const int flowPv ) const
         {
-            const int flowToEbos[ 4 ] = {
+            const int flowToEbos[ 3 ] = {
                                           BlackoilIndices::pressureSwitchIdx,
                                           BlackoilIndices::waterSaturationIdx,
-                                          BlackoilIndices::compositionSwitchIdx,
+                                          BlackoilIndices::compositionSwitchIdx
-                                          BlackoilIndices::solventSaturationIdx
                                         };
+
+            if (flowPv > 2 )
+                return flowPv;
+
             return flowToEbos[ flowPv ];
         }
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const
         {
-            const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
+            const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
+            if (phaseIdx > 2 )
+                return phaseIdx;
             return phaseToComp[ phaseIdx ];
         }
 
@@ -1636,8 +1674,11 @@ namespace Opm {
                     // no need to store refDens for all cells?
                     const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
                     const auto& refDens = intQuants->solventRefDensity();
-                    cellRes[ solventCompIdx ] /= refDens;
+                    cellRes[ contiSolventEqIdx ] /= refDens;
-                    cellRes[ solventCompIdx ] *= cellVolume;
+                    cellRes[ contiSolventEqIdx ] *= cellVolume;
+                }
+                if (has_polymer_) {
+                    cellRes[ contiPolymerEqIdx ] *= cellVolume;
                 }
             }
 
@@ -1670,10 +1711,17 @@ namespace Opm {
                         const auto& refDens = intQuants->solventRefDensity();
                         for( int pvIdx=0; pvIdx < numEq; ++pvIdx )
                         {
-                            (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
+                            (*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
-                            (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
+                            (*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
                         }
                     }
+                    if (has_polymer_) {
+                        for( int pvIdx=0; pvIdx < numEq; ++pvIdx )
+                        {
+                            (*col)[contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
+                        }
+                    }
+
                 }
             }
         }

opm/autodiff/Compat.cpp

@@ -23,7 +23,9 @@
 
 #include <algorithm>
 #include <cassert>
+#include <iostream>
 
+#include <opm/polymer/PolymerBlackoilState.hpp>
 #include <opm/common/data/SimulationDataContainer.hpp>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
@@ -113,6 +115,15 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir,
         sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION );
     }
 
+    if (phases.has_polymer) {
+        if (reservoir.hasCellData( PolymerBlackoilState::CONCENTRATION )) { // compatibility with legacy polymer
+            sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( PolymerBlackoilState::CONCENTRATION ) , data::TargetType::RESTART_SOLUTION );
+        } else {
+            sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::POLYMER ) , data::TargetType::RESTART_SOLUTION );
+        }
+
+    }
+
     return sol;
 }
 

opm/autodiff/FlowMainEbos.hpp

@@ -61,6 +61,7 @@
 namespace Opm
 {
     // The FlowMain class is the ebos based black-oil simulator.
+    template <class TypeTag>
     class FlowMainEbos
     {
         enum FileOutputValue{
@@ -73,7 +74,6 @@ namespace Opm
         };
 
     public:
-        typedef TTAG(EclFlowProblem) TypeTag;
         typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager;
         typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator;
         typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
@@ -82,7 +82,7 @@ namespace Opm
         typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
         typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
 
-        typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator;
+        typedef Opm::SimulatorFullyImplicitBlackoilEbos<TypeTag> Simulator;
         typedef typename Simulator::ReservoirState ReservoirState;
         typedef typename Simulator::OutputWriter OutputWriter;
 
@@ -534,6 +534,27 @@ namespace Opm
             }
 
             initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL"));
+
+            // Get initial polymer concentration from ebos
+            if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
+                auto& cpolymer = state_->getCellData( state_->POLYMER );
+                const int numCells = Opm::UgGridHelpers::numCells(grid);
+                for (int c = 0; c < numCells; ++c) {
+                    cpolymer[c] = ebosProblem().polymerConcentration(c);
+                }
+            }
+            // Get initial solvent saturation from ebos
+            if (GET_PROP_VALUE(TypeTag, EnableSolvent)) {
+                auto& solvent = state_->getCellData( state_->SSOL );
+                auto& sat = state_->saturation();
+                const int np = props.numPhases();
+                const int numCells = Opm::UgGridHelpers::numCells(grid);
+                for (int c = 0; c < numCells; ++c) {
+                    solvent[c] = ebosProblem().solventSaturation(c);
+                    sat[c * np + pu.phase_pos[Water]];
+                }
+            }
+
         }
 
         // Extract messages from parser.
@@ -667,7 +688,7 @@ namespace Opm
         //   fis_solver_
         void setupLinearSolver()
         {
-            typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType;
+            typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
 
             extractParallelGridInformationToISTL(grid(), parallel_information_);
             fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );

opm/autodiff/FlowMainPolymer.hpp

@@ -89,10 +89,6 @@ namespace Opm
             }
         }
 
-
-
-
-
         // Setup linear solver.
         // Writes to:
         //   fis_solver_

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -42,14 +42,12 @@
 
 namespace Opm {
 
-class SimulatorFullyImplicitBlackoilEbos;
-//class StandardWellsDense<FluidSystem>;
 
 /// a simulator for the blackoil model
+template<class TypeTag>
 class SimulatorFullyImplicitBlackoilEbos
 {
 public:
-    typedef typename TTAG(EclFlowProblem) TypeTag;
     typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
     typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@@ -58,11 +56,12 @@ public:
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables)  PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
 
+    typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
 
     typedef WellStateFullyImplicitBlackoilDense WellState;
     typedef BlackoilState ReservoirState;
     typedef BlackoilOutputWriter OutputWriter;
-    typedef BlackoilModelEbos Model;
+    typedef BlackoilModelEbos<TypeTag> Model;
     typedef BlackoilModelParameters ModelParameters;
     typedef NonlinearSolver<Model> Solver;
     typedef StandardWellsDense<TypeTag> WellModel;
@@ -255,7 +254,8 @@ public:
             solver_timer.start();
 
             const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum());
-            WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum());
+            WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_,
+                                 timer.currentStepNum());
 
             auto solver = createSolver(well_model);
 
@@ -413,11 +413,12 @@ public:
     { return ebosSimulator_.gridManager().grid(); }
 
 protected:
-    void handleAdditionalWellInflow(SimulatorTimer& /* timer */,
+    void handleAdditionalWellInflow(SimulatorTimer& timer,
                                     WellsManager& /* wells_manager */,
                                     WellState& /* well_state */,
                                     const Wells* /* wells */)
-    { }
+    {
+    }
 
     std::unique_ptr<Solver> createSolver(WellModel& well_model)
     {
@@ -439,7 +440,8 @@ protected:
                         legacyDepth_,
                         legacyPoreVolume_,
                         rateConverter_.get(),
-                        globalNumCells);
+                        globalNumCells,
+                        grid());
         auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
                                                       model_param_,
                                                       well_model,
@@ -852,6 +854,10 @@ protected:
         }
     }
 
+
+
+
+
     // Data.
     Simulator& ebosSimulator_;
 

opm/autodiff/StandardWellsDense.hpp

@@ -57,6 +57,8 @@
 
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
+#include <math.h>
+
 namespace Opm {
 
 enum WellVariablePositions {
@@ -85,8 +87,11 @@ enum WellVariablePositions {
 
             typedef double Scalar;
             static const int numEq = BlackoilIndices::numEq;
-            static const int numWellEq = numEq; //number of wellEq is the same as numEq in the model
+            static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? 3:numEq; // //numEq; //number of wellEq is only for 3 for polymer
-            static const int solventCompIdx = 3; //TODO get this from ebos
+            static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
+            static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
+            static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
+            static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
             typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
             typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
             typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@@ -94,8 +99,7 @@ enum WellVariablePositions {
             typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
             typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
 
-
+            typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
-
 
             // For the conversion between the surface volume rate and resrevoir voidage rate
             using RateConverterType = RateConverter::
@@ -116,7 +120,8 @@ enum WellVariablePositions {
                       const std::vector<double>& depth_arg,
                       const std::vector<double>& pv_arg,
                       const RateConverterType* rate_converter,
-                      long int global_nc);
+                      long int global_nc,
+                      const auto& grid);
 
 
             /// The number of components in the model.
@@ -294,6 +299,7 @@ enum WellVariablePositions {
             ModelParameters param_;
             bool terminal_output_;
             bool has_solvent_;
+            bool has_polymer_;
             int current_timeIdx_;
 
             PhaseUsage phase_usage_;
@@ -314,6 +320,13 @@ enum WellVariablePositions {
             std::vector<double> well_perforation_densities_;
             std::vector<double> well_perforation_pressure_diffs_;
 
+            std::vector<double> wpolymer_;
+            std::vector<double> wsolvent_;
+
+            std::vector<double> wells_rep_radius_;
+            std::vector<double> wells_perf_length_;
+            std::vector<double> wells_bore_diameter_;
+
             std::vector<EvalWell> wellVariables_;
             std::vector<double> F0_;
 
@@ -397,6 +410,14 @@ enum WellVariablePositions {
                                         const std::vector<double>& initial_potential) const;
 
             double wsolvent(const int well_index) const;
+
+            double wpolymer(const int well_index) const;
+
+            void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
+
+            void computeRepRadiusPerfLength(const auto& grid);
+
+
         };
 
 

opm/autodiff/StandardWellsDense_impl.hpp

@@ -18,6 +18,7 @@ namespace Opm {
        , param_(param)
        , terminal_output_(terminal_output)
        , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
+       , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
        , current_timeIdx_(current_timeIdx)
        , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0)
        , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0)
@@ -47,7 +48,8 @@ namespace Opm {
          const std::vector<double>& depth_arg,
          const std::vector<double>& pv_arg,
          const RateConverterType* rate_converter,
-         long int global_nc)
+         long int global_nc,
+         const auto& grid)
     {
         // has to be set always for the convergence check!
         global_nc_   = global_nc;
@@ -118,6 +120,14 @@ namespace Opm {
         // resize temporary class variables
         Cx_.resize( duneC_.N() );
         invDrw_.resize( invDuneD_.N() );
+
+        if (has_polymer_)
+        {
+            if (PolymerModule::hasPlyshlog()) {
+                computeRepRadiusPerfLength(grid);
+            }
+        }
+
     }
 
 
@@ -197,6 +207,37 @@ namespace Opm {
                 getMobility(ebosSimulator, perf, cell_idx, mob);
                 computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
 
+                if (has_polymer_) {
+                    if (PolymerModule::hasPlyshlog()) {
+                        // compute the well water velocity based on the perforation rates.
+                        double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf];
+                        const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
+                        const auto& scaledDrainageInfo =
+                                materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
+                        const Scalar& Swcr = scaledDrainageInfo.Swcr;
+                        const EvalWell poro = extendEval(intQuants.porosity());
+                        const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
+                        // guard against zero porosity and no water
+                        const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
+                        EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
+                        EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
+
+                        if (PolymerModule::hasShrate()) {
+                            // TODO Use the same conversion as for the reservoar equations.
+                            // Need the "permeability" of the well?
+                            // For now use the same formula as in legacy.
+                            waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf];
+                        }
+                        EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
+                                                                                 intQuants.pvtRegionIndex(),
+                                                                                 waterVelocity);
+
+                        // modify the mobility with the shear factor and recompute the well fluxes.
+                        mob[ Water ] /= shearFactor;
+                        computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
+                    }
+                }
+
                 for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
 
                     // the cq_s entering mass balance equations need to consider the efficiency factors.
@@ -209,17 +250,16 @@ namespace Opm {
                     }
 
                     // subtract sum of phase fluxes in the well equations.
-                    resWell_[w][componentIdx] -= cq_s[componentIdx].value();
+                    resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s[componentIdx].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.
                             ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
-                            duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
+                            duneB_[w][cell_idx][flowToEbosPvIdx(pvIdx)][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
-                            duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
                         }
-                        invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq);
+                        invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq);
                     }
 
                     // add trivial equation for 2p cases (Only support water + oil)
@@ -229,7 +269,7 @@ namespace Opm {
                     }
 
                     // Store the perforation phase flux for later usage.
-                    if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
+                    if (has_solvent_ && componentIdx == solventSaturationIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
                         well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value();
                     } else {
                         well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value();
@@ -455,13 +495,31 @@ namespace Opm {
 
 
 
+    template<typename TypeTag>
+    int
+    StandardWellsDense<TypeTag>::
+    flowToEbosPvIdx( const int flowPv ) const
+    {
+        const int flowToEbos[ 3 ] = {
+            BlackoilIndices::pressureSwitchIdx,
+            BlackoilIndices::waterSaturationIdx,
+            BlackoilIndices::compositionSwitchIdx
+        };
+
+        if (flowPv > 2 )
+            return flowPv;
+
+        return flowToEbos[ flowPv ];
+    }
 
     template<typename TypeTag>
     int
     StandardWellsDense<TypeTag>::
     flowPhaseToEbosCompIdx( const int phaseIdx ) const
     {
-        const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
+        const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
+        if (phaseIdx > 2 )
+            return phaseIdx;
         return phaseToComp[ phaseIdx ];
     }
 
@@ -469,24 +527,6 @@ namespace Opm {
 
 
 
-    template<typename TypeTag>
-    int
-    StandardWellsDense<TypeTag>::
-    flowToEbosPvIdx( const int flowPv ) const
-    {
-        const int flowToEbos[ 4 ] = {
-                                     BlackoilIndices::pressureSwitchIdx,
-                                     BlackoilIndices::waterSaturationIdx,
-                                     BlackoilIndices::compositionSwitchIdx,
-                                     BlackoilIndices::solventSaturationIdx
-                                    };
-        return flowToEbos[ flowPv ];
-    }
-
-
-
-
-
     template<typename TypeTag>
     int
     StandardWellsDense<TypeTag>::
@@ -497,10 +537,6 @@ namespace Opm {
         return flowToEbos[ phaseIdx ];
     }
 
-
-
-
-
     template<typename TypeTag>
     std::vector<double>
     StandardWellsDense<TypeTag>::
@@ -687,8 +723,8 @@ namespace Opm {
     {
         const int nw = wells().number_of_wells;
         // for two-phase numComp < numEq
-        const int numComp = numComponents();
+        //const int numComp = numComponents();
-        for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) {
+        for (int eqIdx = 0; eqIdx < numWellEq;  ++eqIdx) {
             for (int w = 0; w < nw; ++w) {
                 const unsigned int idx = nw * eqIdx + w;
                 assert( idx < wellVariables_.size() );
@@ -697,7 +733,7 @@ namespace Opm {
 
                 eval = 0.0;
                 eval.setValue( xw.wellSolutions()[ idx ] );
-                eval.setDerivative(numWellEq + eqIdx, 1.0);
+                eval.setDerivative(numEq + eqIdx, 1.0);
             }
         }
     }
@@ -765,7 +801,7 @@ namespace Opm {
             b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
         }
         if (has_solvent_) {
-            b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
+            b_perfcells_dense[solventSaturationIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
         }
 
         // Pressure drawdown (also used to determine direction of flow)
@@ -817,7 +853,7 @@ namespace Opm {
             }
 
             if (has_solvent_) {
-                volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx];
+                volumeRatio += cmix_s[solventSaturationIdx] / b_perfcells_dense[solventSaturationIdx];
             }
 
             if (active_[Oil] && active_[Gas]) {
@@ -952,8 +988,9 @@ namespace Opm {
 
         const int nw = wells().number_of_wells;
         const int numComp = numComponents();
-        std::vector<double> res(numComp*nw);
+        std::vector<double> res(numEq*nw, 0.0);
         for( int compIdx = 0; compIdx < numComp; ++compIdx) {
+
             for (int wellIdx = 0; wellIdx < nw; ++wellIdx) {
                 int idx = wellIdx + nw*compIdx;
                 res[idx] = resWell_[ wellIdx ][ compIdx ];
@@ -1006,7 +1043,7 @@ namespace Opm {
                 B += 1 / fs.invB(ebosPhaseIdx).value();
             }
             if (has_solvent_) {
-                auto& B  = B_avg[ solventCompIdx ];
+                auto& B  = B_avg[ solventSaturationIdx ];
                 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
             }
         }
@@ -1204,8 +1241,8 @@ namespace Opm {
 
                 // We use cell values for solvent injector
                 if (has_solvent_) {
-                    b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value();
+                    b_perf[numComp*perf + solventSaturationIdx] = intQuants.solventInverseFormationVolumeFactor().value();
-                    surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity();
+                    surf_dens_perf[numComp*perf + solventSaturationIdx] = intQuants.solventRefDensity();
                 }
             }
         }
@@ -1235,20 +1272,20 @@ namespace Opm {
             // update the second and third well variable (The flux fractions)
             std::vector<double> F(np,0.0);
             if (active_[ Water ]) {
-                const int sign2 = dwells[w][WFrac] > 0 ? 1: -1;
+                const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1;
-                const double dx2_limited = sign2 * std::min(std::abs(dwells[w][WFrac]),dFLimit);
+                const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit);
                 well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited;
             }
 
             if (active_[ Gas ]) {
-                const int sign3 = dwells[w][GFrac] > 0 ? 1: -1;
+                const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1;
-                const double dx3_limited = sign3 * std::min(std::abs(dwells[w][GFrac]),dFLimit);
+                const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit);
                 well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited;
             }
 
             if (has_solvent_) {
-                const int sign4 = dwells[w][SFrac] > 0 ? 1: -1;
+                const int sign4 = dwells[w][flowPhaseToEbosCompIdx(SFrac)] > 0 ? 1: -1;
-                const double dx4_limited = sign4 * std::min(std::abs(dwells[w][SFrac]),dFLimit);
+                const double dx4_limited = sign4 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(SFrac)]),dFLimit);
                 well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited;
             }
 
@@ -1352,7 +1389,7 @@ namespace Opm {
                 case THP: // The BHP and THP both uses the total rate as first well variable.
                 case BHP:
                 {
-                    well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][XvarWell];
+                    well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)];
 
                     switch (wells().type[w]) {
                     case INJECTOR:
@@ -1420,8 +1457,8 @@ namespace Opm {
                 case SURFACE_RATE: // Both rate controls use bhp as first well variable
                 case RESERVOIR_RATE:
                 {
-                    const int sign1 = dwells[w][XvarWell] > 0 ? 1: -1;
+                    const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1;
-                    const double dx1_limited = sign1 * std::min(std::abs(dwells[w][XvarWell]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit);
+                    const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit);
                     well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5);
                     well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w];
 
@@ -1764,7 +1801,7 @@ namespace Opm {
                 perfRates[perf*numComponent + phase] =  xw.perfPhaseRates()[perf*np + phase];
             }
             if(has_solvent_) {
-                perfRates[perf*numComponent + solventCompIdx] =  xw.perfRateSolvent()[perf];
+                perfRates[perf*numComponent + solventSaturationIdx] =  xw.perfRateSolvent()[perf];
             }
         }
         well_perforation_densities_ =
@@ -2143,7 +2180,7 @@ namespace Opm {
         if (wells().type[wellIdx] == INJECTOR) {
             if (has_solvent_ ) {
                 double comp_frac = 0.0;
-                if (compIdx == solventCompIdx) { // solvent
+                if (has_solvent_ && compIdx == solventSaturationIdx) { // solvent
                     comp_frac = wells().comp_frac[np*wellIdx + pu.phase_pos[ Gas ]] * wsolvent(wellIdx);
                 } else if (compIdx == pu.phase_pos[ Gas ]) {
                     comp_frac = wells().comp_frac[np*wellIdx + compIdx] * (1.0 - wsolvent(wellIdx));
@@ -2209,7 +2246,7 @@ namespace Opm {
                 EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(wellIdx, phase_under_control);
                 if (has_solvent_ && phase_under_control == Gas) {
                     // for GRAT controlled wells solvent is included in the target
-                    wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventCompIdx);
+                    wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventSaturationIdx);
                 }
 
                 if (compIdx == phase_under_control) {
@@ -2274,7 +2311,7 @@ namespace Opm {
             return wellVariables_[GFrac * nw + wellIdx];
         }
 
-        if (compIdx == solventCompIdx) {
+        if (has_solvent_ && compIdx == solventSaturationIdx) {
             return wellVariables_[SFrac * nw + wellIdx];
         }
 
@@ -2305,7 +2342,7 @@ namespace Opm {
         const WellControls* wc = wells().ctrls[wellIdx];
         if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
 
-            if (has_solvent_ && compIdx == solventCompIdx) {
+            if (has_solvent_ && compIdx == solventSaturationIdx) {
                 return wellVolumeFraction(wellIdx, compIdx);
             }
             const double* distr = well_controls_get_current_distr(wc);
@@ -2710,10 +2747,10 @@ namespace Opm {
                 xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
             }
             if (active_[ Gas ]) {
-                xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
+                xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (xw.wellRates()[np*well_index + Gas] - xw.solventWellRate(well_index)) / tot_well_rate ;
             }
             if (has_solvent_) {
-                xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * wsolvent(well_index) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
+                xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * xw.solventWellRate(well_index) / tot_well_rate ;
             }
         } else {
             const WellType& well_type = wells().type[well_index];
@@ -3027,6 +3064,171 @@ namespace Opm {
     }
 
 
+    template<typename TypeTag>
+    double
+    StandardWellsDense<TypeTag>::
+    wpolymer(const int well_index) const {
+
+        if (!has_polymer_) {
+            return 0.0;
+        }
+
+        // loop over all wells until we find the well with the matching name
+        for (const auto&  well : wells_ecl_) {
+            if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) {
+                continue;
+            }
+
+            WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_);
+            WellPolymerProperties polymer = well->getPolymerProperties(current_timeIdx_);
+            if (injection.injectorType == WellInjector::WATER) {
+
+                double polymerFraction = polymer.m_polymerConcentration;
+
+                // Look until we find the correct well
+                if (well->name() == wells().name[well_index]) {
+                    return polymerFraction;
+                }
+            }
+        }
+        // we didn't find it return 0;
+        assert(false);
+        return 0.0;
+    }
+
+
+    template<typename TypeTag>
+    void
+    StandardWellsDense<TypeTag>::
+    setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const
+    {
+        if (global_cell) {
+            for (int i = 0; i < number_of_cells; ++i) {
+                cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
+            }
+        }
+        else {
+            for (int i = 0; i < number_of_cells; ++i) {
+                cartesian_to_compressed.insert(std::make_pair(i, i));
+            }
+        }
+
+    }
+
+    template<typename TypeTag>
+    void
+    StandardWellsDense<TypeTag>::
+    computeRepRadiusPerfLength(const auto& grid)
+    {
+
+        // TODO, the function does not work for parallel running
+        // to be fixed later.
+        int number_of_cells = Opm::UgGridHelpers::numCells(grid);
+        const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
+        const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
+        auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
+        auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
+
+        if (wells_ecl_.size() == 0) {
+            OPM_MESSAGE("No wells specified in Schedule section, "
+                        "initializing no wells");
+            return;
+        }
+
+        const int nw = wells().number_of_wells;
+        const int nperf = wells().well_connpos[nw];
+
+        const size_t timeStep = current_timeIdx_;
+
+        wells_rep_radius_.clear();
+        wells_perf_length_.clear();
+        wells_bore_diameter_.clear();
+
+        wells_rep_radius_.reserve(nperf);
+        wells_perf_length_.reserve(nperf);
+        wells_bore_diameter_.reserve(nperf);
+
+        std::map<int,int> cartesian_to_compressed;
+
+        setupCompressedToCartesian(global_cell, number_of_cells,
+                                    cartesian_to_compressed);
+
+        int well_index = 0;
+
+        for (auto wellIter= wells_ecl_.begin(); wellIter != wells_ecl_.end(); ++wellIter) {
+             const auto* well = (*wellIter);
+
+             if (well->getStatus(timeStep) == WellCommon::SHUT) {
+                 continue;
+             }
+             {   // COMPDAT handling
+                 const auto& completionSet = well->getCompletions(timeStep);
+                 for (size_t c=0; c<completionSet.size(); c++) {
+                     const auto& completion = completionSet.get(c);
+                     if (completion.getState() == WellCompletion::OPEN) {
+                         int i = completion.getI();
+                         int j = completion.getJ();
+                         int k = completion.getK();
+
+                         const int* cpgdim = cart_dims;
+                         int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
+                         std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
+                         if (cgit == cartesian_to_compressed.end()) {
+                             OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
+                                       << k << " not found in grid (well = " << well->name() << ')');
+                         }
+                         int cell = cgit->second;
+
+                         {
+                             double radius = 0.5*completion.getDiameter();
+                             if (radius <= 0.0) {
+                                 radius = 0.5*unit::feet;
+                                 OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
+                             }
+
+                             const std::array<double, 3> cubical =
+                             WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
+
+                             WellCompletion::DirectionEnum direction = completion.getDirection();
+
+                             double re; // area equivalent radius of the grid block
+                             double perf_length; // the length of the well perforation
+
+                             switch (direction) {
+                                 case Opm::WellCompletion::DirectionEnum::X:
+                                     re = std::sqrt(cubical[1] * cubical[2] / M_PI);
+                                     perf_length = cubical[0];
+                                     break;
+                                 case Opm::WellCompletion::DirectionEnum::Y:
+                                     re = std::sqrt(cubical[0] * cubical[2] / M_PI);
+                                     perf_length = cubical[1];
+                                     break;
+                                 case Opm::WellCompletion::DirectionEnum::Z:
+                                     re = std::sqrt(cubical[0] * cubical[1] / M_PI);
+                                     perf_length = cubical[2];
+                                     break;
+                                 default:
+                                     OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
+                             }
+
+                             double repR = std::sqrt(re * radius);
+                             wells_rep_radius_.push_back(repR);
+                             wells_perf_length_.push_back(perf_length);
+                             wells_bore_diameter_.push_back(2. * radius);
+                         }
+                     } else {
+                         if (completion.getState() != WellCompletion::SHUT) {
+                             OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
+                         }
+                     }
+
+                 }
+            }
+            well_index++;
+        }
+    }
+
+
 
 
 } // namespace Opm