substitute solveWellEq with one inner solve for newly opened wells

ebos/ebos.hh

@@ -207,7 +207,6 @@ public:
         EWOMS_HIDE_PARAM(TypeTag, MaxInnerIterWells);
         EWOMS_HIDE_PARAM(TypeTag, MaxSinglePrecisionDays);
         EWOMS_HIDE_PARAM(TypeTag, MaxStrictIter);
-        EWOMS_HIDE_PARAM(TypeTag, SolveWelleqInitially);
         EWOMS_HIDE_PARAM(TypeTag, UpdateEquationsScaling);
         EWOMS_HIDE_PARAM(TypeTag, UseUpdateStabilization);
         EWOMS_HIDE_PARAM(TypeTag, MatrixAddWellContributions);

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -84,10 +84,6 @@ struct MaxStrictIter {
     using type = UndefinedProperty;
 };
 template<class TypeTag, class MyTypeTag>
-struct SolveWelleqInitially {
-    using type = UndefinedProperty;
-};
-template<class TypeTag, class MyTypeTag>
 struct UpdateEquationsScaling {
     using type = UndefinedProperty;
 };
@@ -213,10 +209,6 @@ struct MaxStrictIter<TypeTag, TTag::FlowModelParameters> {
     static constexpr int value = 0;
 };
 template<class TypeTag>
-struct SolveWelleqInitially<TypeTag, TTag::FlowModelParameters> {
-    static constexpr bool value = true;
-};
-template<class TypeTag>
 struct UpdateEquationsScaling<TypeTag, TTag::FlowModelParameters> {
     static constexpr bool value = false;
 };
@@ -361,9 +353,6 @@ namespace Opm
         /// Maximum number of Newton iterations before we give up on the CNV convergence criterion
         int max_strict_iter_;
 
-        /// Solve well equation initially
-        bool solve_welleq_initially_;
-
         /// Update scaling factors for mass balance equations
         bool update_equations_scaling_;
 
@@ -409,7 +398,6 @@ namespace Opm
             max_inner_iter_wells_ = EWOMS_GET_PARAM(TypeTag, int, MaxInnerIterWells);
             maxSinglePrecisionTimeStep_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays) *24*60*60;
             max_strict_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxStrictIter);
-            solve_welleq_initially_ = EWOMS_GET_PARAM(TypeTag, bool, SolveWelleqInitially);
             update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling);
             use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization);
             matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
@@ -444,7 +432,6 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, RegularizationFactorMsw, "Regularization factor for ms wells");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays, "Maximum time step size where single precision floating point arithmetic can be used solving for the linear systems of equations");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIter, "Maximum number of Newton iterations before relaxed tolerances are used for the CNV convergence criterion");
-            EWOMS_REGISTER_PARAM(TypeTag, bool, SolveWelleqInitially, "Fully solve the well equations before each iteration of the reservoir model");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UpdateEquationsScaling, "Update scaling factors for mass balance equations during the run");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method");
             EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices");

opm/simulators/timestepping/SimulatorReport.cpp

@@ -156,17 +156,7 @@ namespace Opm
 
         }
 
-        int n = total_well_iterations + (failureReport ? failureReport->total_well_iterations : 0);
+        int n = total_linearizations + (failureReport ? failureReport->total_linearizations : 0);
-
-        os << fmt::format("Overall Well Iterations:   {:7}", n);
-        if (failureReport) {
-          os << fmt::format("    (Failed: {:3}; {:2.1f}%)",
-                            failureReport->total_well_iterations,
-                            100.0*failureReport->total_well_iterations/n);
-        }
-        os << std::endl;
-
-        n = total_linearizations + (failureReport ? failureReport->total_linearizations : 0);
         os << fmt::format("Overall Linearizations:    {:7}", n);
         if (failureReport) {
           os << fmt::format("    (Failed: {:3}; {:2.1f}%)",

opm/simulators/wells/BlackoilWellModel.hpp

@@ -415,8 +415,6 @@ namespace Opm {
             /// at the beginning of the time step and no derivatives are included in these quantities
             void calculateExplicitQuantities(Opm::DeferredLogger& deferred_logger) const;
 
-            SimulatorReportSingle solveWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
-
             void initPrimaryVariablesEvaluation() const;
 
             // The number of components in the model.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -406,6 +406,17 @@ namespace Opm {
         //compute well guideRates
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         WellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, well_state_, comm, guideRate_.get());
+        updateAndCommunicateGroupData();
+        // Compute initial well solution for new wells
+        for (auto& well : well_container_) {
+            const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE;
+            const bool event = report_step_starts_ && schedule().hasWellGroupEvent(well->name(), effective_events_mask, reportStepIdx);
+            if (event) {
+                well->calculateExplicitQuantities(ebosSimulator_, well_state_, local_deferredLogger);
+                well->solveWellToInitialize(ebosSimulator_, well_state_, local_deferredLogger);
+            }
+        }
+
         logAndCheckForExceptionsAndThrow(local_deferredLogger,
             exception_thrown, "beginTimeStep() failed.", terminal_output_);
 
@@ -951,23 +962,6 @@ namespace Opm {
 
             std::vector< Scalar > B_avg(numComponents(), Scalar() );
             computeAverageFormationFactor(B_avg);
-
-            if (param_.solve_welleq_initially_ && iterationIdx == 0) {
-                // solve the well equations as a pre-processing step
-                last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
-
-
-                if (initial_step_) {
-                    // update the explicit quantities to get the initial fluid distribution in the well correct.
-                    calculateExplicitQuantities(local_deferredLogger);
-                    prepareTimeStep(local_deferredLogger);
-                    last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
-                    initial_step_ = false;
-                }
-                // TODO: should we update the explicit related here again, or even prepareTimeStep().
-                // basically, this is a more updated state from the solveWellEq based on fixed
-                // reservoir state, will tihs be a better place to inialize the explict information?
-            }
             gliftDebug("assemble() : running assembleWellEq()..", local_deferredLogger);
             well_state_.enableGliftOptimization();
             assembleWellEq(B_avg, dt, local_deferredLogger);
@@ -1163,84 +1157,6 @@ namespace Opm {
 
 
 
-    template<typename TypeTag>
-    SimulatorReportSingle
-    BlackoilWellModel<TypeTag>::
-    solveWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger)
-    {
-        WellState well_state0 = well_state_;
-
-        const int max_iter = param_.max_welleq_iter_;
-
-        int it  = 0;
-        bool converged;
-        int exception_thrown = 0;
-        do {
-            try {
-                assembleWellEq(B_avg, dt, deferred_logger);
-            } catch (std::exception& e) {
-                exception_thrown = 1;
-            }
-            // We need to check on all processes, as getWellConvergence() below communicates on all processes.
-            logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
-
-            const auto report = getWellConvergence(B_avg);
-            converged = report.converged();
-
-            if (converged) {
-                break;
-            }
-
-            try {
-                if( localWellsActive() )
-                {
-                    for (auto& well : well_container_) {
-                        well->solveEqAndUpdateWellState(well_state_, deferred_logger);
-                    }
-                }
-                // updateWellControls uses communication
-                // Therefore the following is executed if there
-                // are active wells anywhere in the global domain.
-                if( wellsActive() )
-                {
-                    updateWellControls(deferred_logger, /*don't switch group controls*/false);
-                    initPrimaryVariablesEvaluation();
-                }
-            } catch (std::exception& e) {
-                exception_thrown = 1;
-            }
-
-            logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
-            ++it;
-        } while (it < max_iter);
-
-        try {
-            if (converged) {
-                if (terminal_output_) {
-                    deferred_logger.debug("Well equation solution gets converged with " + std::to_string(it) + " iterations");
-                }
-            } else {
-                if (terminal_output_) {
-                    deferred_logger.debug("Well equation solution failed in getting converged with " + std::to_string(it) + " iterations");
-                }
-                well_state_ = well_state0;
-                updatePrimaryVariables(deferred_logger);
-            }
-        } catch (std::exception& e) {
-            exception_thrown = 1;
-        }
-
-        logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
-
-        SimulatorReportSingle report;
-        report.converged = converged;
-        report.total_well_iterations = it;
-        return report;
-    }
-
-
-
-
 
     template<typename TypeTag>
     ConvergenceReport

opm/simulators/wells/WellInterface.hpp

@@ -320,6 +320,11 @@ namespace Opm
 
         void setDynamicThpLimit(const double thp_limit);
 
+        void solveWellToInitialize(const Simulator& ebosSimulator,
+                                   WellState& well_state,
+                                   Opm::DeferredLogger& deferred_logger);
+
+
     protected:
 
         // to indicate a invalid completion

opm/simulators/wells/WellInterface_impl.hpp

@@ -1385,6 +1385,27 @@ namespace Opm
         }
     }
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    solveWellToInitialize(const Simulator& ebosSimulator,
+                          WellState& well_state,
+                          Opm::DeferredLogger& deferred_logger)
+    {
+        // keep a copy of the original well state
+        const WellState well_state0 = well_state;
+        const double dt = ebosSimulator.timeStepSize();
+        const bool converged = iterateWellEquations(ebosSimulator, B_avg_, dt, well_state, deferred_logger);
+        if (converged) {
+            deferred_logger.debug("Compute initial well solution for well " + name() +  ". Converged");
+        } else {
+            const int max_iter = param_.max_welleq_iter_;
+            deferred_logger.debug("Compute initial well solution for well " +name() + ". Failed to converge in "
+                                  + std::to_string(max_iter) + " iterations");
+            well_state = well_state0;
+        }
+    }
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const