trying to add an inner well iteration to speed up convergence.

opm/autodiff/BlackoilModelParameters.cpp

@@ -53,6 +53,7 @@ namespace Opm
         tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_);
         max_welleq_iter_ = param.getDefault("max_welleq_iter", max_welleq_iter_);
         max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_);
+        use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_);
         maxSinglePrecisionTimeStep_ = unit::convert::from(
                 param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day );
         max_strict_iter_ = param.getDefault("max_strict_iter",8);
@@ -81,7 +82,8 @@ namespace Opm
         tolerance_well_control_ = 1.0e-7;
         tolerance_pressure_ms_wells_ = 100.0;
         max_welleq_iter_ = 15;
-        max_pressure_change_ms_wells_ = 100000.; // 1.0 bar
+        max_pressure_change_ms_wells_ = 200000.; // 2.0 bar
+        use_inner_iterations_ms_wells_ = true;
         maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day );
         solve_welleq_initially_ = true;
         update_equations_scaling_ = false;

opm/autodiff/BlackoilModelParameters.hpp

@@ -56,6 +56,9 @@ namespace Opm
         /// Maximum pressure change over an iteratio for ms wells
         double max_pressure_change_ms_wells_;
 
+        /// Whether to use inner iterations for ms wells
+        bool use_inner_iterations_ms_wells_;
+
         /// Maximum iteration number of the well equation solution
         int max_welleq_iter_;
 

opm/autodiff/BlackoilWellModel.hpp

@@ -212,7 +212,7 @@ namespace Opm {
             void computeRepRadiusPerfLength(const Grid& grid);
 
 
-            void computeAverageFormationFactor(Simulator& ebosSimulator,
+            void computeAverageFormationFactor(const Simulator& ebosSimulator,
                                                std::vector<double>& B_avg) const;
 
             void applyVREPGroupControl(WellState& well_state) const;

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -210,7 +210,7 @@ namespace Opm {
                    bool only_wells) const
     {
         for (int w = 0; w < number_of_wells_; ++w) {
-            well_container_[w]->assembleWellEq(ebosSimulator, dt, well_state, only_wells);
+            well_container_[w]->assembleWellEq(ebosSimulator, param_, dt, well_state, only_wells);
         }
     }
 
@@ -929,7 +929,7 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    computeAverageFormationFactor(Simulator& ebosSimulator,
+    computeAverageFormationFactor(const Simulator& ebosSimulator,
                                   std::vector<double>& B_avg) const
     {
         const int np = numPhases();

opm/autodiff/MultisegmentWell.hpp

@@ -105,6 +105,7 @@ namespace Opm
         virtual void initPrimaryVariablesEvaluation() const;
 
         virtual void assembleWellEq(Simulator& ebosSimulator,
+                                    const ModelParameters& param,
                                     const double dt,
                                     WellState& well_state,
                                     bool only_wells);
@@ -115,7 +116,7 @@ namespace Opm
                                                WellState& well_state) const;
 
         /// check whether the well equations get converged for this well
-        virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
+        virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
                                                      const std::vector<double>& B_avg,
                                                      const ModelParameters& param) const;
 
@@ -327,6 +328,17 @@ namespace Opm
         bool frictionalPressureLossConsidered() const;
 
         bool accelerationalPressureLossConsidered() const;
+
+        // TODO: try to make ebosSimulator const, as it should be
+        void iterateWellEquations(Simulator& ebosSimulator,
+                                  const ModelParameters& param,
+                                  const double dt,
+                                  WellState& well_state);
+
+        void assembleWellEqWithoutIteration(Simulator& ebosSimulator,
+                                            const double dt,
+                                            WellState& well_state,
+                                            bool only_wells);
     };
 
 }

opm/autodiff/MultisegmentWell_impl.hpp

@@ -220,120 +220,18 @@ namespace Opm
     void
     MultisegmentWell<TypeTag>::
     assembleWellEq(Simulator& ebosSimulator,
+                   const ModelParameters& param,
                    const double dt,
                    WellState& well_state,
                    bool only_wells)
     {
-        // clear all entries
+
-        if (!only_wells) {
+        const bool use_inner_iterations = param.use_inner_iterations_ms_wells_;
-            duneB_ = 0.0;
+        if (use_inner_iterations) {
-            duneC_ = 0.0;
+            iterateWellEquations(ebosSimulator, param, dt, well_state);
         }
 
-        duneD_ = 0.0;
+        assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
-        resWell_ = 0.0;
-
-        // for the black oil cases, there will be four equations,
-        // the first three of them are the mass balance equations, the last one is the pressure equations.
-        //
-        // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
-
-        auto& ebosJac = ebosSimulator.model().linearizer().matrix();
-        auto& ebosResid = ebosSimulator.model().linearizer().residual();
-
-        const bool allow_cf = getAllowCrossFlow();
-
-        const int nseg = numberOfSegments();
-        const int num_comp = numComponents();
-
-        // TODO: finding better place to put it
-        computeSegmentFluidProperties(ebosSimulator);
-
-        for (int seg = 0; seg < nseg; ++seg) {
-            // calculating the accumulation term // TODO: without considering the efficiencty factor for now
-            // volume of the segment
-            {
-                const double volume = segmentSet()[seg].volume();
-                // for each component
-                for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
-                    EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx])
-                                               + getSegmentRate(seg, comp_idx);
-
-                    resWell_[seg][comp_idx] += accumulation_term.value();
-                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-                        duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
-                    }
-                }
-            }
-
-            // considering the contributions from the inlet segments
-            {
-                for (const int inlet : segment_inlets_[seg]) {
-                    for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
-                        const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx);
-                        resWell_[seg][comp_idx] -= inlet_rate.value();
-                        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-                            duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq);
-                        }
-                    }
-                }
-            }
-
-            // calculating the perforation rate for each perforation that belongs to this segment
-            const EvalWell seg_pressure = getSegmentPressure(seg);
-            for (const int perf : segment_perforations_[seg]) {
-                const int cell_idx = well_cells_[perf];
-                const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
-                std::vector<EvalWell> mob(num_comp, 0.0);
-                getMobility(ebosSimulator, perf, mob);
-                std::vector<EvalWell> cq_s(num_comp, 0.0);
-                computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s);
-
-                for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
-                    // the cq_s entering mass balance equations need to consider the efficiency factors.
-                    const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
-
-                    if (!only_wells) {
-                        // subtract sum of component fluxes in the reservoir equation.
-                        // need to consider the efficiency factor
-                        // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input
-                        // is a component index :D
-                        ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
-                    }
-
-                    // subtract sum of phase fluxes in the well equations.
-                    resWell_[seg][comp_idx] -= cq_s_effective.value();
-
-                    // assemble the jacobians
-                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-                        if (!only_wells) {
-                            // also need to consider the efficiency factor when manipulating the jacobians.
-                            duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
-                        }
-                        // the index name for the D should be eq_idx / pv_idx
-                        duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
-                    }
-
-                    for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
-                        if (!only_wells) {
-                            // also need to consider the efficiency factor when manipulating the jacobians.
-                            ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx);
-                            duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
-                        }
-                    }
-                }
-                // TODO: we should save the perforation pressure and preforation rates?
-                // we do not use it in the simulation for now, while we might need them if
-                // we handle the pressure in SEG mode.
-            }
-
-            // the fourth dequation, the pressure drop equation
-            if (seg == 0) { // top segment, pressure equation is the control equation
-                assembleControlEq();
-            } else {
-                assemblePressureEq(seg);
-            }
-        }
     }
 
 
@@ -491,7 +389,7 @@ namespace Opm
     template <typename TypeTag>
     typename MultisegmentWell<TypeTag>::ConvergenceReport
     MultisegmentWell<TypeTag>::
-    getWellConvergence(Simulator& ebosSimulator,
+    getWellConvergence(const Simulator& /* ebosSimulator */,
                        const std::vector<double>& B_avg,
                        const ModelParameters& param) const
     {
@@ -823,6 +721,10 @@ namespace Opm
                     const BlackoilModelParameters& param,
                     WellState& well_state) const
     {
+        const bool use_inner_iterations = param.use_inner_iterations_ms_wells_;
+
+        const double relaxation_factor = use_inner_iterations ? 0.4 : 1.0;
+
         // I guess the following can also be applied to the segmnet pressure
         // maybe better to give it a different name
         const double dBHPLimit = param.dbhp_max_rel_;
@@ -833,13 +735,13 @@ namespace Opm
         for (int seg = 0; seg < numberOfSegments(); ++seg) {
             if (active()[ Water ]) {
                 const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
-                const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), dFLimit);
+                const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit);
                 primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
             }
 
             if (active()[ Gas ]) {
                 const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
-                const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), dFLimit);
+                const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit);
                 primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
             }
 
@@ -856,7 +758,7 @@ namespace Opm
 
             // update the total rate // TODO: should we have a limitation of the total rate change?
             {
-                primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - dwells[seg][GTotal];
+                primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];
             }
 
             // TODO: not handling solvent related for now
@@ -1860,4 +1762,177 @@ namespace Opm
         return (segmentSet().compPressureDrop() == WellSegment::HFA);
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    iterateWellEquations(Simulator& ebosSimulator,
+                         const ModelParameters& param,
+                         const double dt,
+                         WellState& well_state)
+    {
+        std::cout << " beginning iterateWellEquations " << std::endl;
+        // basically, it only iterate through the equations.
+        // we update the primary variables
+        // if converged, we can update the well_state.
+        // the function updateWellState() should have a flag to show
+        // if we will update the well state.
+        // assembleWellEq(
+        const int max_iter_number = 3;
+        int it = 0;
+        for (; it < max_iter_number; ++it) {
+
+            assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true);
+
+            const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_);
+
+            // TODO: use these small values for now, not intend to reach the convergence
+            // in this stage, but, should we?
+            // If we want to use the real one, we need to find a way to get them.
+            const std::vector<double> B {0.1, 0.1, 0.001};
+
+            const ConvergenceReport report = getWellConvergence(ebosSimulator, B, param);
+
+            if (report.converged) {
+                std::cout << " converged in iterateWellEquations " << std::endl;
+                break;
+            }
+
+            updateWellState(dx_well, param, well_state);
+
+            initPrimaryVariablesEvaluation();
+        }
+
+        if (it >= max_iter_number) {
+            std::cout << " the iterateWellEquations did not converged " << std::endl;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    assembleWellEqWithoutIteration(Simulator& ebosSimulator,
+                                   const double dt,
+                                   WellState& well_state,
+                                   bool only_wells)
+    {
+        // calculate the fluid properties needed.
+        computeSegmentFluidProperties(ebosSimulator);
+
+        // clear all entries
+        if (!only_wells) {
+            duneB_ = 0.0;
+            duneC_ = 0.0;
+        }
+
+        duneD_ = 0.0;
+        resWell_ = 0.0;
+
+        // for the black oil cases, there will be four equations,
+        // the first three of them are the mass balance equations, the last one is the pressure equations.
+        //
+        // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
+
+        auto& ebosJac = ebosSimulator.model().linearizer().matrix();
+        auto& ebosResid = ebosSimulator.model().linearizer().residual();
+
+        const bool allow_cf = getAllowCrossFlow();
+
+        const int nseg = numberOfSegments();
+        const int num_comp = numComponents();
+
+        for (int seg = 0; seg < nseg; ++seg) {
+            // calculating the accumulation term // TODO: without considering the efficiencty factor for now
+            // volume of the segment
+            {
+                const double volume = segmentSet()[seg].volume();
+                // for each component
+                for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
+                    EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx])
+                                               + getSegmentRate(seg, comp_idx);
+
+                    resWell_[seg][comp_idx] += accumulation_term.value();
+                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+                        duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
+                    }
+                }
+            }
+
+            // considering the contributions from the inlet segments
+            {
+                for (const int inlet : segment_inlets_[seg]) {
+                    for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
+                        const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx);
+                        resWell_[seg][comp_idx] -= inlet_rate.value();
+                        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+                            duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq);
+                        }
+                    }
+                }
+            }
+
+            // calculating the perforation rate for each perforation that belongs to this segment
+            const EvalWell seg_pressure = getSegmentPressure(seg);
+            for (const int perf : segment_perforations_[seg]) {
+                const int cell_idx = well_cells_[perf];
+                const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
+                std::vector<EvalWell> mob(num_comp, 0.0);
+                getMobility(ebosSimulator, perf, mob);
+                std::vector<EvalWell> cq_s(num_comp, 0.0);
+                computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s);
+
+                for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
+                    // the cq_s entering mass balance equations need to consider the efficiency factors.
+                    const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;
+
+                    if (!only_wells) {
+                        // subtract sum of component fluxes in the reservoir equation.
+                        // need to consider the efficiency factor
+                        // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input
+                        // is a component index :D
+                        ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value();
+                    }
+
+                    // subtract sum of phase fluxes in the well equations.
+                    resWell_[seg][comp_idx] -= cq_s_effective.value();
+
+                    // assemble the jacobians
+                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+                        if (!only_wells) {
+                            // also need to consider the efficiency factor when manipulating the jacobians.
+                            duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix
+                        }
+                        // the index name for the D should be eq_idx / pv_idx
+                        duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq);
+                    }
+
+                    for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
+                        if (!only_wells) {
+                            // also need to consider the efficiency factor when manipulating the jacobians.
+                            ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(comp_idx)][pv_idx] -= cq_s_effective.derivative(pv_idx);
+                            duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx);
+                        }
+                    }
+                }
+                // TODO: we should save the perforation pressure and preforation rates?
+                // we do not use it in the simulation for now, while we might need them if
+                // we handle the pressure in SEG mode.
+            }
+
+            // the fourth dequation, the pressure drop equation
+            if (seg == 0) { // top segment, pressure equation is the control equation
+                assembleControlEq();
+            } else {
+                assemblePressureEq(seg);
+            }
+        }
+    }
+
 }

opm/autodiff/StandardWell.hpp

@@ -111,6 +111,7 @@ namespace Opm
         virtual void initPrimaryVariablesEvaluation() const;
 
         virtual void assembleWellEq(Simulator& ebosSimulator,
+                                    const ModelParameters& param,
                                     const double dt,
                                     WellState& well_state,
                                     bool only_wells);
@@ -121,7 +122,7 @@ namespace Opm
                                                WellState& xw) const;
 
         /// check whether the well equations get converged for this well
-        virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
+        virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
                                                      const std::vector<double>& B_avg,
                                                      const ModelParameters& param) const;
 

opm/autodiff/StandardWell_impl.hpp

@@ -514,6 +514,7 @@ namespace Opm
     void
     StandardWell<TypeTag>::
     assembleWellEq(Simulator& ebosSimulator,
+                   const ModelParameters& /* param */,
                    const double dt,
                    WellState& well_state,
                    bool only_wells)
@@ -1364,7 +1365,7 @@ namespace Opm
     template<typename TypeTag>
     typename StandardWell<TypeTag>::ConvergenceReport
     StandardWell<TypeTag>::
-    getWellConvergence(Simulator& ebosSimulator,
+    getWellConvergence(const Simulator& /* ebosSimulator */,
                        const std::vector<double>& B_avg,
                        const ModelParameters& param) const
     {

opm/autodiff/WellInterface.hpp

@@ -144,7 +144,7 @@ namespace Opm
             }
         };
 
-        virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
+        virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
                                                      const std::vector<double>& B_avg,
                                                      const ModelParameters& param) const = 0;
 
@@ -152,6 +152,7 @@ namespace Opm
                                                WellState& well_state) = 0;
 
         virtual void assembleWellEq(Simulator& ebosSimulator,
+                                    const ModelParameters& param,
                                     const double dt,
                                     WellState& well_state,
                                     bool only_wells) = 0;