Solve well equation for the first iteration

ebos/ebos.hh

@@ -207,6 +207,7 @@ 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);

flow/flow_ebos_gasoil.cpp

@@ -38,6 +38,10 @@ struct EclFlowGasOilProblem {
     using InheritsFrom = std::tuple<EclFlowProblem>;
 };
 }
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowGasOilProblem> {
+    static constexpr bool value = true;
+};
 
 //! The indices required by the model
 template<class TypeTag>

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -84,6 +84,10 @@ 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;
 };
@@ -209,6 +213,10 @@ 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;
 };
@@ -353,6 +361,9 @@ 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_;
 
@@ -398,6 +409,7 @@ 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);
@@ -432,6 +444,7 @@ 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/wells/BlackoilWellModel_impl.hpp

@@ -409,11 +409,12 @@ namespace Opm {
         updateAndCommunicateGroupData();
         // Compute initial well solution for new wells
         for (auto& well : well_container_) {
-            const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE;
+            const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+                                                 + ScheduleEvents::NEW_WELL;
             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);
+                well->solveWellEquation(ebosSimulator_, well_state_, local_deferredLogger);
             }
         }
 
@@ -962,6 +963,14 @@ namespace Opm {
 
             std::vector< Scalar > B_avg(numComponents(), Scalar() );
             computeAverageFormationFactor(B_avg);
+
+            if (param_.solve_welleq_initially_ && iterationIdx == 0) {
+                for (auto& well : well_container_) {
+                    well->solveWellEquation(ebosSimulator_, well_state_, local_deferredLogger);
+                }
+                updateWellControls(local_deferredLogger, /* check group controls */ true);
+            }
+
             gliftDebug("assemble() : running assembleWellEq()..", local_deferredLogger);
             well_state_.enableGliftOptimization();
             assembleWellEq(B_avg, dt, local_deferredLogger);

opm/simulators/wells/WellInterface.hpp

@@ -320,9 +320,9 @@ namespace Opm
 
         void setDynamicThpLimit(const double thp_limit);
 
-        void solveWellToInitialize(const Simulator& ebosSimulator,
+        void solveWellEquation(const Simulator& ebosSimulator,
-                                   WellState& well_state,
+                               WellState& well_state,
-                                   Opm::DeferredLogger& deferred_logger);
+                               Opm::DeferredLogger& deferred_logger);
 
 
     protected:

opm/simulators/wells/WellInterface_impl.hpp

@@ -1388,9 +1388,9 @@ namespace Opm
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
-    solveWellToInitialize(const Simulator& ebosSimulator,
+    solveWellEquation(const Simulator& ebosSimulator,
-                          WellState& well_state,
+                      WellState& well_state,
-                          Opm::DeferredLogger& deferred_logger)
+                      Opm::DeferredLogger& deferred_logger)
     {
         // keep a copy of the original well state
         const WellState well_state0 = well_state;