Merge pull request #1970 from atgeirr/handle-grup-unititialized

[WIP] Handle GRUP-caused uninitialized (-1) well control index

CMakeLists_files.cmake

@@ -32,7 +32,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/core/wells/WellCollection.cpp
   opm/core/wells/WellsGroup.cpp
   opm/core/wells/WellsManager.cpp
-  opm/core/wells/well_controls.c
+  opm/core/wells/well_controls.cpp
   opm/core/wells/wells.c
   opm/simulators/flow/MissingFeatures.cpp
   opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp

opm/core/wells/well_controls.cpp

@@ -42,6 +42,8 @@
 #include <stdlib.h>
 #include <stdio.h>
 
+#include <stdexcept>
+
 /**
  * Controls for a single well.
  * Each control specifies a well rate or bottom-hole pressure. Only
@@ -142,7 +144,7 @@ well_controls_create(void)
 {
     struct WellControls *ctrl;
 
-    ctrl = malloc(1 * sizeof *ctrl);
+    ctrl = static_cast<WellControls*>(malloc(1 * sizeof *ctrl));
 
     if (ctrl != NULL) {
         /* Initialise empty control set; the well is created open. */
@@ -176,11 +178,11 @@ well_controls_reserve(int nctrl, struct WellControls *ctrl)
     distr  = realloc(ctrl->distr , nctrl * ctrl->number_of_phases * sizeof *ctrl->distr );
 
     int ok = 0;
-    if (type   != NULL) { ctrl->type   = type  ; ok++; }
+    if (type   != NULL) { ctrl->type   = static_cast<WellControlType*>(type)  ; ok++; }
-    if (target != NULL) { ctrl->target = target; ok++; }
+    if (target != NULL) { ctrl->target = static_cast<double*>(target); ok++; }
-    if (alq    != NULL) { ctrl->alq    = alq;    ok++; }
+    if (alq    != NULL) { ctrl->alq    = static_cast<double*>(alq   ); ok++; }
-    if (vfp    != NULL) { ctrl->vfp    = vfp;    ok++; }
+    if (vfp    != NULL) { ctrl->vfp    = static_cast<int*>(vfp);       ok++; }
-    if (distr  != NULL) { ctrl->distr  = distr ; ok++; }
+    if (distr  != NULL) { ctrl->distr  = static_cast<double*>(distr) ; ok++; }
 
     if (ok == 5) {
         for (int c = ctrl->cpty; c < nctrl; c++) {
@@ -204,18 +206,18 @@ struct WellControls *
 well_controls_clone(const struct WellControls *ctrl)
 /* ---------------------------------------------------------------------- */
 {
-    struct WellControls* new = well_controls_create();
+    struct WellControls* new_ctrls = well_controls_create();
 
-    if (new != NULL) {
+    if (new_ctrls != NULL) {
         /* Assign appropriate number of phases */
-        well_controls_assert_number_of_phases(new, ctrl->number_of_phases);
+        well_controls_assert_number_of_phases(new_ctrls, ctrl->number_of_phases);
 
         int n  = well_controls_get_num(ctrl);
-        int ok = well_controls_reserve(n, new);
+        int ok = well_controls_reserve(n, new_ctrls);
 
         if (! ok) {
-            well_controls_destroy(new);
+            well_controls_destroy(new_ctrls);
-            new = NULL;
+            new_ctrls= NULL;
         }
         else {
             int i;
@@ -226,32 +228,32 @@ well_controls_clone(const struct WellControls *ctrl)
                 double alq = well_controls_iget_alq   (ctrl, i);
                 int vfp = well_controls_iget_vfp   (ctrl, i);
 
-                ok = well_controls_add_new(type, target, alq, vfp, distr, new);
+                ok = well_controls_add_new(type, target, alq, vfp, distr, new_ctrls);
             }
 
             if (i < n) {
                 assert (!ok);
-                well_controls_destroy(new);
+                well_controls_destroy(new_ctrls);
 
-                new = NULL;
+                new_ctrls = NULL;
             }
             else {
                 i = well_controls_get_current(ctrl);
-                well_controls_set_current(new, i);
+                well_controls_set_current(new_ctrls, i);
 
                 if (well_controls_well_is_open(ctrl)) {
-                    well_controls_open_well(new);
+                    well_controls_open_well(new_ctrls);
                 }
                 else {
-                    well_controls_stop_well(new);
+                    well_controls_stop_well(new_ctrls);
                 }
             }
         }
     }
 
-    assert (well_controls_equal(ctrl, new, true));
+    assert (well_controls_equal(ctrl, new_ctrls, true));
 
-    return new;
+    return new_ctrls;
 }
 
 
@@ -295,6 +297,8 @@ well_controls_iget_type(const struct WellControls * ctrl, int control_index) {
 
 enum WellControlType 
 well_controls_get_current_type(const struct WellControls * ctrl) {
+    if (ctrl->current < 0)
+        throw std::logic_error("Tried to use invalid current control < 0");
     return well_controls_iget_type( ctrl , ctrl->current);
 }
 
@@ -312,6 +316,8 @@ well_controls_iget_target(const struct WellControls * ctrl, int control_index) {
 
 double
 well_controls_get_current_target(const struct WellControls * ctrl) {
+    if (ctrl->current < 0)
+        throw std::logic_error("Tried to use invalid current control < 0");
     return ctrl->target[ctrl->current];
 }
 
@@ -350,6 +356,8 @@ well_controls_iget_distr(const struct WellControls * ctrl, int control_index) {
 
 const double *
 well_controls_get_current_distr(const struct WellControls * ctrl) {
+    if (ctrl->current < 0)
+        throw std::logic_error("Tried to use invalid current control < 0");
     return well_controls_iget_distr( ctrl , ctrl->current );
 }
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -19,6 +19,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/WellInjectionProperties.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
 namespace Opm
@@ -2506,12 +2507,8 @@ namespace Opm
             }
         }
 
-
+        if (std::abs(total_well_rate) > 0.) {
-        const WellControls* wc = well_controls_;
-        const double* distr = well_controls_get_current_distr(wc);
             const auto pu = phaseUsage();
-
-        if(std::abs(total_well_rate) > 0.) {
             if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                 primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
             }
@@ -2523,9 +2520,10 @@ namespace Opm
             }
         } else { // total_well_rate == 0
             if (well_type_ == INJECTOR) {
+                auto phase = well_ecl_.getInjectionProperties().injectorType;
                 // only single phase injection handled
                 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    if (distr[Water] > 0.0) {
+                    if (phase == WellInjector::TypeEnum::WATER) {
                         primary_variables_[WFrac] = 1.0;
                     } else {
                         primary_variables_[WFrac] = 0.0;
@@ -2533,7 +2531,7 @@ namespace Opm
                 }
 
                 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    if (distr[pu.phase_pos[Gas]] > 0.0) {
+                    if (phase == WellInjector::TypeEnum::GAS) {
                         primary_variables_[GFrac] = 1.0 - wsolvent();
                         if (has_solvent) {
                             primary_variables_[SFrac] = wsolvent();

opm/simulators/wells/WellInterface_impl.hpp

@@ -1205,9 +1205,8 @@ namespace Opm
     WellInterface<TypeTag>::scalingFactor(const int phaseIdx) const
     {
         const WellControls* wc = well_controls_;
-        const double* distr = well_controls_get_current_distr(wc);
 
-        if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
+        if (well_controls_get_current(wc) != -1 && well_controls_get_current_type(wc) == RESERVOIR_RATE) {
             if (has_solvent && phaseIdx == contiSolventEqIdx ) {
                 typedef Ewoms::BlackOilSolventModule<TypeTag> SolventModule;
                 double coeff = 0;
@@ -1215,6 +1214,7 @@ namespace Opm
                 return coeff;
             }
             // TODO: use the rateConverter here as well.
+            const double* distr = well_controls_get_current_distr(wc);
             return distr[phaseIdx];
         }
         const auto& pu = phaseUsage();

opm/simulators/wells/WellState.hpp

@@ -109,6 +109,19 @@ namespace Opm
                             const int first_cell = wells->well_cells[wells->well_connpos[w]];
                             bhp_[w] = cellPressures[first_cell];
                         }
+                    } else if (well_controls_get_current(ctrl) == -1) {
+                        // Well under group control.
+                        // 1. Rates: assign zero well rates.
+                        for (int p = 0; p < np; ++p) {
+                            wellrates_[np*w + p] = 0.0;
+                        }
+                        // 2. Bhp: initialize bhp to be a
+                        //    little above or below (depending on if
+                        //    the well is an injector or producer)
+                        //    pressure in first perforation cell.
+                        const int first_cell = wells->well_cells[wells->well_connpos[w]];
+                        const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
+                        bhp_[w] = safety_factor*cellPressures[first_cell];
                     } else {
                         // Open well:
                         // 1. Rates: initialize well rates to match controls