adding flow_injecitivity executable target

CMakeLists.txt

@@ -129,7 +129,8 @@ opm_add_test(flow
     flow/flow_ebos_polymer.cpp
     flow/flow_ebos_solvent.cpp
     flow/flow_ebos_energy.cpp
-    flow/flow_ebos_oilwater_polymer.cpp)
+    flow/flow_ebos_oilwater_polymer.cpp
+    flow/flow_ebos_oilwater_polymer_injectivity.cpp)
 install(TARGETS flow DESTINATION bin)
 
 add_test(NAME flow__version

flow/flow_ebos_oilwater_polymer_injectivity.cpp

@@ -0,0 +1,98 @@
+/*
+  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/>.
+*/
+#include "config.h"
+
+// Define making clear that the simulator supports AMG
+#define FLOW_SUPPORT_AMG 1
+
+#include <flow/flow_ebos_oilwater_polymer_injectivity.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <ewoms/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/autodiff/FlowMainEbos.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/misc/mpimanager.hh>
+#else
+#include <dune/common/parallel/mpihelper.hh>
+#endif
+
+namespace Ewoms {
+namespace Properties {
+NEW_TYPE_TAG(EclFlowOilWaterPolymerProblem, INHERITS_FROM(EclFlowProblem));
+SET_BOOL_PROP(EclFlowOilWaterPolymerProblem, EnablePolymer, true);
+SET_BOOL_PROP(EclFlowOilWaterPolymerProblem, EnablePolymerMW, true);
+//! The indices required by the model
+// For this case, there will be two primary variables introduced for the polymer
+// polymer concentration and polymer molecular weight
+// TODO: probaby it can be better to refer to the implementation of flow_ebos_oilwater, or other two phase
+// simulators. Not sure why they make it more complicated.
+/* SET_TYPE_PROP(EclFlowOilWaterPolymerProblem, Indices,
+    Ewoms::BlackOilTwoPhaseIndices< 0,
+              2,
+              0,
+              0,
+              2>); */
+SET_PROP(EclFlowOilWaterPolymerProblem, Indices)
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    typedef TTAG(EclFlowProblem) BaseTypeTag;
+    typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem;
+
+public:
+    typedef Ewoms::BlackOilTwoPhaseIndices<0,
+                                           2,
+                                           0,
+                                           /*PVOffset=*/0,
+                                           /*disabledCompIdx=*/FluidSystem::gasCompIdx> type;
+};
+}}
+
+namespace Opm {
+/* void flowEbosOilWaterPolymerInjectivitySetDeck(Deck& deck, EclipseState& eclState, Schedule& schedule, SummaryConfig& summaryConfig)
+{
+    typedef TTAG(EclFlowOilWaterPolymerProblem) TypeTag;
+    typedef GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
+
+    Vanguard::setExternalDeck(&deck, &eclState, &schedule, &summaryConfig);
+} */
+
+// ----------------- Main program -----------------
+int flowEbosOilWaterPolymerInjectivityMain(int argc, char** argv)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    Opm::resetLocale();
+
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc, argv);
+#else
+    Dune::MPIHelper::instance(argc, argv);
+#endif
+
+    Opm::FlowMainEbos<TTAG(EclFlowOilWaterPolymerProblem)> mainfunc;
+
+    return mainfunc.execute(argc, argv);
+}
+
+}

flow/flow_ebos_oilwater_polymer_injectivity.hpp

@@ -0,0 +1,31 @@
+/*
+  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 FLOW_EBOS_OILWATER_POLYMER_INJECTIVITY_HPP
+#define FLOW_EBOS_OILWATER_POLYMER_INJECTIVITY_HPP
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
+
+
+namespace Opm {
+  // void flowEbosOilWaterPolymerInjectivitySetDeck(Deck& deck, EclipseState& eclState, Schedule& schedule, SummaryConfig& summary_config);
+  int flowEbosOilWaterPolymerInjectivityMain(int argc, char** argv);
+}
+
+#endif // FLOW_EBOS_OILWATER_POLYMER_INJECTIVITY_HPP

flow/flow_polymer_injectivity.cpp

@@ -0,0 +1,25 @@
+/*
+  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/>.
+*/
+#include "config.h"
+
+#include <flow/flow_ebos_oilwater_polymer_injectivity.hpp>
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    return Opm::flowEbosOilWaterPolymerInjectivityMain(argc, argv);
+}

opm/autodiff/BlackoilModelEbos.hpp

@@ -120,8 +120,10 @@ namespace Opm {
         static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
         static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
         static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
+        static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
         static const int solventSaturationIdx = Indices::solventSaturationIdx;
         static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
+        static const int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
         static const int temperatureIdx = Indices::temperatureIdx;
 
         typedef Dune::FieldVector<Scalar, numEq >        VectorBlockType;
@@ -157,6 +159,7 @@ namespace Opm {
         , has_vapoil_(FluidSystem::enableVaporizedOil())
         , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
         , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
+        , has_polymermw_(GET_PROP_VALUE(TypeTag, EnablePolymerMW))
         , has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy))
         , param_( param )
         , well_model_ (well_model)
@@ -767,6 +770,19 @@ namespace Opm {
                     R_sum[ contiPolymerEqIdx ] += R2;
                     maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
                 }
+
+                if (has_polymermw_) {
+                    assert(has_polymer_);
+
+                    B_avg[contiPolymerMWEqIdx] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
+                    // the residual of the polymer molecular equatinon is scaled down  by a 100, since molecular weight
+                    // can be much bigger than 1, and this equation shares the same tolerance with other mass balance equations
+                    // TODO: there should be a more general way to determine the scaling-down coefficient
+                    const auto R2 = ebosResid[cell_idx][contiPolymerMWEqIdx] / 100.;
+                    R_sum[contiPolymerMWEqIdx] += R2;
+                    maxCoeff[contiPolymerMWEqIdx] = std::max( maxCoeff[contiPolymerMWEqIdx], std::abs( R2 ) / pvValue );
+                }
+
                 if (has_energy_ ) {
                     B_avg[ contiEnergyEqIdx ] += 1.0;
                     const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
@@ -833,6 +849,10 @@ namespace Opm {
                 if (has_polymer_) {
                     compNames[polymerConcentrationIdx] = "Polymer";
                 }
+                if (has_polymermw_) {
+                    assert(has_polymer_);
+                    compNames[polymerMoleWeightIdx] = "PolymerMW";
+                }
                 if (has_energy_) {
                     compNames[temperatureIdx] = "Energy";
                 }
@@ -993,6 +1013,7 @@ namespace Opm {
         const bool has_vapoil_;
         const bool has_solvent_;
         const bool has_polymer_;
+        const bool has_polymermw_;
         const bool has_energy_;
 
         ModelParameters                 param_;

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -146,13 +146,6 @@ namespace Opm
                             perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
                         }
                         perfPress()[perf] = cellPressures[wells->well_cells[perf]];
-
-                        // here we did not consider the case that we close some perforation during the running
-                        // and also, wells can be shut and re-opened
-                        // it should only be updated when necessary, while WellState might not have the facility to do it
-                        perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ];
-                        perf_skin_pressure_[ perf ] = prevState->perfSkinPressure()[ oldPerf_idx ];
-                        perf_water_velocity_[ perf ] = prevState->perfWaterVelocity()[ oldPerf_idx ];
                     }
                 }
             }
@@ -243,6 +236,23 @@ namespace Opm
                                 }
                             }
                         }
+
+                        // polymer injectivity related
+                        // here we did not consider the case that we close some perforation during the running
+                        // and also, wells can be shut and re-opened
+                        if (pu.has_polymermw) {
+                            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 )
+                                {
+                                    perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ];
+                                    perf_skin_pressure_[ perf ] = prevState->perfSkinPressure()[ oldPerf_idx ];
+                                    perf_water_velocity_[ perf ] = prevState->perfWaterVelocity()[ oldPerf_idx ];
+                                }
+                            }
+                        }
                     }
 
 

opm/core/props/phaseUsageFromDeck.hpp

@@ -108,7 +108,7 @@ namespace Opm
             pu.phase_pos[BlackoilPhases::Energy] = -1;
 
         // Add polymer molecular weight related
-        pu.has_polymermw = phase.active(Phase::POLYMW);
+        pu.has_polymermw = phases.active(Phase::POLYMW);
         if (pu.has_polymermw) {
             if (!pu.has_polymer) {
                 OPM_THROW(std::runtime_error, "pu.has_polymermw is true while pu.has_polymer is false");