calculating the representative radius and perf length

for all well perforations, to be used in shear-thinning calculation.

The calculation is approximated.

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -80,6 +80,8 @@ namespace Opm {
                              const bool                              has_vapoil,
                              const bool                              has_polymer,
                              const bool                              has_plyshlog,
+                             std::vector<double>&                    wells_rep_radius,
+                             std::vector<double>&                    wells_perf_length,
                              const bool                              terminal_output);
 
         /// Called once before each time step.
@@ -128,6 +130,11 @@ namespace Opm {
         const int  poly_pos_;
         V cmax_;
 
+        // representative radius and perforation length of well perforations
+        // to be used in shear-thinning computation.
+        std::vector<double> wells_rep_radius_;
+        std::vector<double> wells_perf_length_;
+
         // Need to declare Base members we want to use here.
         using Base::grid_;
         using Base::fluid_;
@@ -252,6 +259,9 @@ namespace Opm {
                              std::vector<double>& maxNormWell,
                              int nc,
                              int nw) const;
+
+
+
     };
 
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -86,13 +86,17 @@ namespace Opm {
                                                      const bool                              has_vapoil,
                                                      const bool                              has_polymer,
                                                      const bool                              has_plyshlog,
+                                                     std::vector<double>&                    wells_rep_radius,
+                                                     std::vector<double>&                    wells_perf_length,
                                                      const bool                              terminal_output)
         : Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
                has_disgas, has_vapoil, terminal_output),
           polymer_props_ad_(polymer_props_ad),
           has_polymer_(has_polymer),
           has_plyshlog_(has_plyshlog),
-          poly_pos_(detail::polymerPos(fluid.phaseUsage()))
+          poly_pos_(detail::polymerPos(fluid.phaseUsage())),
+          wells_rep_radius_(wells_rep_radius),
+          wells_perf_length_(wells_perf_length)
     {
         if (has_polymer_) {
             if (!active_[Water]) {
@@ -576,10 +580,6 @@ namespace Opm {
         return converged;
     }
 
-
-
-
-
     template <class Grid>
     ADB
     BlackoilPolymerModel<Grid>::computeMc(const SolutionState& state) const
@@ -587,8 +587,6 @@ namespace Opm {
         return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
     }
 
-
-
 } // namespace Opm
 
 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp

@@ -124,17 +124,37 @@ namespace Opm
 
         std::unique_ptr<Solver> createSolver(const Wells* wells);
 
+
         void handleAdditionalWellInflow(SimulatorTimer& timer,
                                         WellsManager& wells_manager,
                                         typename BaseType::WellState& well_state,
                                         const Wells* wells);
 
+
     private:
         const PolymerPropsAd& polymer_props_;
         bool has_polymer_;
-        /// flag for PLYSHLOG keywords
+        // flag for PLYSHLOG keyword
         bool has_plyshlog_;
         DeckConstPtr deck_;
+
+        std::vector<double> wells_rep_radius_;
+        std::vector<double> wells_perf_length_;
+
+        // generate the mapping from Cartesian grid cells to global compressed cells,
+        // copied from opm-core, to be used in function computeRepRadiusPerfLength()
+        static void
+        setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed);
+
+        //  calculate the representative radius and length for for well peforations
+        //  it will be used in the shear-thinning calcluation only.
+        void
+        computeRepRadiusPerfLength(const Opm::EclipseStateConstPtr eclipseState,
+                                   const size_t                    timeStep,
+                                   const GridT&                    grid,
+                                   std::vector<double>&            wells_rep_radius,
+                                   std::vector<double>&            wells_perf_length);
+
     };
 
 } // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp

@@ -68,6 +68,7 @@ namespace Opm
         ModelParams modelParams( BaseType::param_ );
         typedef NewtonSolver<Model> Solver;
 
+
         auto model = std::unique_ptr<Model>(new Model(modelParams,
                                                       BaseType::grid_,
                                                       BaseType::props_,
@@ -80,6 +81,8 @@ namespace Opm
                                                       BaseType::has_vapoil_,
                                                       has_polymer_,
                                                       has_plyshlog_,
+                                                      wells_rep_radius_,
+                                                      wells_perf_length_,
                                                       BaseType::terminal_output_));
 
         if (!BaseType::threshold_pressures_by_face_.empty()) {
@@ -91,6 +94,9 @@ namespace Opm
         return std::unique_ptr<Solver>(new Solver(solverParams, std::move(model)));
     }
 
+
+
+
     template <class GridT>
     void SimulatorFullyImplicitBlackoilPolymer<GridT>::
     handleAdditionalWellInflow(SimulatorTimer& timer,
@@ -115,6 +121,134 @@ namespace Opm
                                             timer.simulationTimeElapsed() + timer.currentStepLength(),
                                             polymer_inflow_c);
         well_state.polymerInflow() = polymer_inflow_c;
+
+        computeRepRadiusPerfLength(BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_);
+    }
+
+
+    template <class GridT>
+    void SimulatorFullyImplicitBlackoilPolymer<GridT>::
+    setupCompressedToCartesian(const int* global_cell, int number_of_cells,
+                               std::map<int,int>& cartesian_to_compressed )
+    {
+        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 <class GridT>
+    void SimulatorFullyImplicitBlackoilPolymer<GridT>::
+    computeRepRadiusPerfLength(const Opm::EclipseStateConstPtr eclipseState,
+                               const size_t                    timeStep,
+                               const GridT&                    grid,
+                               std::vector<double>&            wells_rep_radius,
+                               std::vector<double>&            wells_perf_length)
+    {
+
+        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 (eclipseState->getSchedule()->numWells() == 0) {
+            OPM_MESSAGE("No wells specified in Schedule section, "
+                        "initializing no wells");
+            return;
+        }
+
+        wells_rep_radius.clear();
+        wells_perf_length.clear();
+
+        std::map<int,int> cartesian_to_compressed;
+
+        setupCompressedToCartesian(global_cell, number_of_cells,
+                                   cartesian_to_compressed);
+
+        ScheduleConstPtr          schedule = eclipseState->getSchedule();
+        std::vector<WellConstPtr> wells    = schedule->getWells(timeStep);
+
+        int well_index = 0;
+
+        for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
+             WellConstPtr well = (*wellIter);
+
+             if (well->getStatus(timeStep) == WellCommon::SHUT) {
+                 continue;
+             }
+             {   // COMPDAT handling
+                 CompletionSetConstPtr completionSet = well->getCompletions(timeStep);
+                 for (size_t c=0; c<completionSet->size(); c++) {
+                     CompletionConstPtr 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);
+                         }
+                     } else {
+                         if (completion->getState() != WellCompletion::SHUT) {
+                             OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion->getState() ) << " not handled");
+                         }
+                     }
+
+                 }
+            }
+            well_index++;
+        }
     }
 
 } // namespace Opm