Merge pull request #3000 from totto82/substituteSolveWellEq

substitute solveWellEq with inner iterations and initial solve

compareECLFiles.cmake

@@ -950,7 +950,7 @@ if(MPI_FOUND)
                                        FILENAME SPE1CASE1_BRINE
                                        SIMULATOR flow
                                        ABS_TOL ${abs_tol_parallel}
-                                       REL_TOL ${rel_tol_parallel}
+                                       REL_TOL ${coarse_rel_tol_parallel}
                                        TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
 
   add_test_compare_parallel_simulation(CASENAME fetkovich_2d

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,30 @@ namespace Opm {
         //compute well guideRates
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         WellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, well_state_, comm, guideRate_.get());
+        try {
+            updateAndCommunicateGroupData();
+            // Compute initial well solution for new wells
+            for (auto& well : well_container_) {
+                const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+                                                    + ScheduleEvents::NEW_WELL;
+
+                const auto& events = schedule()[reportStepIdx].wellgroup_events();
+                const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
+                if (event) {
+                    try {
+                        well->calculateExplicitQuantities(ebosSimulator_, well_state_, local_deferredLogger);
+                        well->solveWellEquation(ebosSimulator_, well_state_, local_deferredLogger);
+                    } catch (std::exception& e) {
+                        const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
+                        local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
+                    }
+                }
+            }
+        } catch (std::exception& e) {
+            const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
+            local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
+        }
+
         logAndCheckForExceptionsAndThrow(local_deferredLogger,
             exception_thrown, "beginTimeStep() failed.", terminal_output_);
 
@@ -953,21 +977,12 @@ namespace Opm {
             computeAverageFormationFactor(B_avg);
 
             if (param_.solve_welleq_initially_ && iterationIdx == 0) {
-                // solve the well equations as a pre-processing step
+                for (auto& well : well_container_) {
-                last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
+                    well->solveWellEquation(ebosSimulator_, well_state_, 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().
+                updateWellControls(local_deferredLogger, /* check group controls */ false);
-                // 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 +1178,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

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

opm/simulators/wells/WellInterface_impl.hpp

@@ -1402,6 +1402,30 @@ namespace Opm
         }
     }
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    solveWellEquation(const Simulator& ebosSimulator,
+                      WellState& well_state,
+                      Opm::DeferredLogger& deferred_logger)
+    {
+        if (!this->isOperable())
+            return;
+
+        // 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