Merge pull request #972 from GitPaean/fixing_polymer_omega

fixing the calculation when polymer not fully mixed

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -227,6 +227,10 @@ namespace Opm {
         computeAccum(const SolutionState& state,
                      const int            aix  );
 
+        void
+        computeInjectionMobility(const SolutionState& state,
+                                 std::vector<ADB>& mob_perfcells);
+
         void
         assembleMassBalanceEq(const SolutionState& state);
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -266,7 +266,6 @@ namespace Opm {
                         concentration.begin() ,
                         max_concentration.begin() ,
                         [](double c_max , double c) { return std::max( c_max , c ); });
-
     }
 
 
@@ -463,14 +462,13 @@ namespace Opm {
             if (has_polymer_) {
                 const ADB tr_mult = transMult(state.pressure);
                 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
-                const ADB mc = computeMc(state);
                 const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
                 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
                 // Reduce mobility of water phase by relperm reduction and effective viscosity increase.
                 sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
                 // Compute polymer mobility.
                 const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
-                sd_.rq[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc;
+                sd_.rq[poly_pos_].mob = tr_mult * krw_eff * state.concentration * inv_poly_eff_visc;
                 sd_.rq[poly_pos_].b = sd_.rq[actph].b;
                 sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
                 UpwindSelector<double> upwind(grid_, ops_, sd_.rq[poly_pos_].dh.value());
@@ -541,6 +539,11 @@ namespace Opm {
         std::vector<ADB> mob_perfcells;
         std::vector<ADB> b_perfcells;
         wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
+
+        // updating the the injection mobility related to polymer injection when necessary
+        // only the mobility of water phase is updated
+        computeInjectionMobility(state, mob_perfcells);
+
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
             Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
@@ -566,9 +569,9 @@ namespace Opm {
             V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
             ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v);
             mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
+            wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
         }
 
-        wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
         wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
         wellModel().addWellFluxEq(cq_s, state, residual_);
         addWellContributionToMassBalanceEq(cq_s, state, well_state);
@@ -622,7 +625,6 @@ namespace Opm {
         UpwindSelector<double> upwind(grid_, ops_, dh.value());
 
         const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
-        const ADB mc = computeMc(state);
         ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
                                                          cmax,
                                                          sd_.rq[phase].kr);
@@ -752,6 +754,15 @@ namespace Opm {
             if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
                 visc_mult_wells[i] = 1.;
             }
+            if (selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
+                if (xw.polymerInflow()[well_cells[i]] == 0.) {
+                    visc_mult_wells[i] = 1.;
+                } else {
+                    // TODO: not tested for this assumption yet
+                    const double c_perf = state.concentration.value()[well_cells[i]];
+                    visc_mult_wells[i] = polymer_props_ad_.viscMult(c_perf);
+                }
+            }
         }
 
         const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
@@ -778,6 +789,71 @@ namespace Opm {
         return;
     }
 
+
+
+
+
+    template<class Grid>
+    void
+    BlackoilPolymerModel<Grid>::
+    computeInjectionMobility(const SolutionState& state,
+                             std::vector<ADB>& mob_perfcells)
+    {
+        // Calculating the mobility for the polymer injection peforations
+        if (has_polymer_ && wellModel().localWellsActive()) {
+            const std::vector<int> well_cells = wellModel().wellOps().well_cells;
+            const int nperf = well_cells.size();
+
+            // Calculating the drawdown to decide the injection perforation
+            const ADB& p_perfcells = subset(state.pressure, well_cells);
+            const V& cdp = wellModel().wellPerforationPressureDiffs();
+            const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
+            // Pressure drawdown (also used to determine direction of flow)
+            const ADB drawdown =  p_perfcells - perfpressure;
+
+            // Polymer concentration in the perforations
+            const ADB c_perfcells = subset(state.concentration, well_cells);
+
+            // Distinguishing the injection perforation from other perforation
+            // The value is the location in the well_cell array, not the global index
+            std::vector<int> polymer_inj_cells;
+            std::vector<int> other_well_cells;
+
+            polymer_inj_cells.reserve(nperf);
+            other_well_cells.reserve(nperf);
+
+            for (int c = 0; c < nperf; ++c) {
+                // TODO: more tests need to be done for this criterion
+                if (drawdown.value()[c] < 0.0 && c_perfcells.value()[c] > 0.0) {
+                    polymer_inj_cells.push_back(c);
+                } else {
+                    other_well_cells.push_back(c);
+                }
+            }
+
+            // there is some polymer injection process going
+            if ( !polymer_inj_cells.empty() ) {
+                // the mobility need to be recalculated for the polymer injection cells
+                const int water_pos = fluid_.phaseUsage().phase_pos[Water];
+                const ADB mu_perfcells = subset(sd_.rq[water_pos].mu, well_cells);
+
+                const ADB c_poly_inj_cells = subset(c_perfcells, polymer_inj_cells);
+                const ADB mu_poly_inj_cells = subset(mu_perfcells, polymer_inj_cells); // water viscosity
+
+                const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(c_poly_inj_cells, mu_poly_inj_cells.value());
+                const ADB fully_mixing_visc = polymer_props_ad_.viscMult(c_poly_inj_cells) * mu_poly_inj_cells;
+
+                // the original mobility for the polymer injection well cells
+                ADB mob_polymer_inj = subset(mob_perfcells[water_pos], polymer_inj_cells);
+                const ADB mob_others = subset(mob_perfcells[water_pos], other_well_cells);
+
+                mob_polymer_inj = mob_polymer_inj / inv_wat_eff_visc / fully_mixing_visc;
+                mob_perfcells[water_pos] = superset(mob_polymer_inj, polymer_inj_cells, nperf) +
+                                           superset(mob_others, other_well_cells, nperf);
+            }
+        }
+    }
+
 } // namespace Opm
 
 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED

opm/polymer/fullyimplicit/PolymerPropsAd.cpp

@@ -130,6 +130,29 @@ namespace Opm {
     }
 
 
+    ADB
+    PolymerPropsAd::viscMult(const ADB& c) const
+    {
+        const int nc = c.size();
+        V visc_mult(nc);
+        V dvisc_mult(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            double im = 0, dim = 0;
+            im = polymer_props_.viscMultWithDer(c.value()(i), &dim);
+            visc_mult(i) = im;
+            dvisc_mult(i) = dim;
+        }
+
+        ADB::M dim_diag(dvisc_mult.matrix().asDiagonal());
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dim_diag * c.derivative()[block];
+        }
+        return ADB::function(std::move(visc_mult), std::move(jacs));
+    }
+
 
 
     PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props)

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -74,6 +74,9 @@ namespace Opm {
         /// \param[in] c        Array of n polymer concentraion values.
         /// \return             Array of n viscosity multiplier from PLVISC table.
 
+
+        ADB viscMult(const ADB& c) const;
+
         /// Constructor wrapping a polymer props.
         PolymerPropsAd(const PolymerProperties& polymer_props);