added test for consitency

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1819,6 +1819,20 @@ namespace Opm
                                    const GroupState<Scalar>& group_state,
                                    DeferredLogger& deferred_logger)
     {
+        //auto value = this->primary_variables_;
+        auto value = this->getPrimaryVars();
+        this->updatePrimaryVariables(simulator, well_state, deferred_logger);
+        auto new_value = this->getPrimaryVars();
+        assert(new_value.size() == value.size());
+        for(int i=0; i < new_value.size(); ++i){
+            assert(value[i] == new_value[i]);
+        }
+        // for(int i=0; i<this->numberOfSegments(); ++i) {
+        //     for(int j =0; j <Indices::numEq; ++j) {
+        //         assert(value[i][j] == this->primary_variables_[i][j]);
+        //     }
+        // }
+
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // update the upwinding segments

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -117,6 +117,104 @@ resize(const int numWellEq)
     numWellEq_ = numWellEq;
 }
 
+template<class FluidSystem, class Indices>
+bool StandardWellPrimaryVariables<FluidSystem,Indices>::
+consistent(const WellState<Scalar>& well_state) const
+{
+    static constexpr int Water = BlackoilPhases::Aqua;
+    static constexpr int Oil = BlackoilPhases::Liquid;
+    static constexpr int Gas = BlackoilPhases::Vapour;
+
+    const int well_index = well_.indexOfWell();
+    const int np = well_.numPhases();
+    const auto& pu = well_.phaseUsage();
+    const auto& ws = well_state.well(well_index);
+    // the weighted total well rate
+    Scalar total_well_rate = 0.0;
+    for (int p = 0; p < np; ++p) {
+        total_well_rate += well_.scalingFactor(p) * ws.surface_rates[p];
+    }
+
+    // Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
+    // under surface condition is used here
+    auto approx = [](auto a,auto b){ return std::abs(a-b) < 1e-10; };
+    bool ok = true;
+    if (well_.isInjector()) {
+        switch (well_.wellEcl().injectorType()) {
+        case InjectorType::WATER:
+            ok = ok && approx(value_[WQTotal],ws.surface_rates[pu.phase_pos[Water]]);
+            break;
+        case InjectorType::GAS:
+            ok = ok && approx(value_[WQTotal],ws.surface_rates[pu.phase_pos[Gas]]);
+            break;
+        case InjectorType::OIL:
+            ok = ok && approx(value_[WQTotal], ws.surface_rates[pu.phase_pos[Oil]]);
+            break;
+        case InjectorType::MULTI:
+            // Not supported.
+            ok = false;
+            assert(false);
+            break;
+        }
+    } else {
+        ok = ok && approx(value_[WQTotal], total_well_rate);     
+    }
+    if (std::abs(total_well_rate) > 0.) {
+        if constexpr (has_wfrac_variable) {
+            auto tmp = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate; 
+            ok = ok && approx(value_[WFrac], tmp);
+        }
+        if constexpr (has_gfrac_variable) {
+            auto tmp =  well_.scalingFactor(pu.phase_pos[Gas]) *
+                            (ws.surface_rates[pu.phase_pos[Gas]] -
+                             (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate;
+            ok = ok && approx(value_[GFrac], tmp) ;
+        }
+        if constexpr (Indices::enableSolvent) {
+            //value_[SFrac] = well_.scalingFactor(Indices::contiSolventEqIdx) * ws.sum_solvent_rates() / total_well_rate ;
+        }
+
+    } else { // total_well_rate == 0
+        if (well_.isInjector()) {
+            // only single phase injection handled
+            if constexpr (has_wfrac_variable) {
+                // if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                //     auto phase = well_.wellEcl().getInjectionProperties().injectorType;
+                //     if (phase == InjectorType::WATER) {
+                //         value_[WFrac] = 1.0;
+                //     } else {
+                //         value_[WFrac] = 0.0;
+                //     }
+                // }
+            }
+            if constexpr (has_gfrac_variable) {
+                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                    auto phase = well_.wellEcl().getInjectionProperties().injectorType;
+                    if (phase == InjectorType::GAS) {
+                        ok = ok && approx(value_[GFrac], (1.0 - well_.rsRvInj()));
+                      
+                    } else {
+                        ok = ok && approx(value_[GFrac], 0.0);
+                    }
+                }
+            }
+
+            // TODO: it is possible to leave injector as a oil well,
+            // when F_w and F_g both equals to zero, not sure under what kind of circumstance
+            // this will happen.
+        } else if (well_.isProducer()) { // producers
+       
+        } else {
+           assert(false);
+        }
+    }
+
+    // BHP
+    ok = ok && approx(value_[Bhp], ws.bhp);
+    return ok;
+}
+
+
 template<class FluidSystem, class Indices>
 void StandardWellPrimaryVariables<FluidSystem,Indices>::
 update(const WellState<Scalar>& well_state,

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -101,6 +101,9 @@ public:
     //! \brief Returns number of well equations.
     int numWellEq() const { return numWellEq_; }
 
+    bool consistent(const WellState<Scalar>& well_state) const;
+
+
     //! \brief Copy values from well state.
     void update(const WellState<Scalar>& well_state,
                 const bool stop_or_zero_rate_target,

opm/simulators/wells/StandardWell_impl.hpp

@@ -356,7 +356,8 @@ namespace Opm
                                    const GroupState<Scalar>& group_state,
                                    DeferredLogger& deferred_logger)
     {
-        // TODO: only_wells should be put back to save some computation
+        assert(this->primary_variables_.consistent(well_state));
+         // TODO: only_wells should be put back to save some computation
         // for example, the matrices B C does not need to update if only_wells
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 

opm/simulators/wells/WellInterface_impl.hpp

@@ -845,11 +845,8 @@ namespace Opm
                                 DeferredLogger& deferred_logger)
     {
         OPM_TIMEFUNCTION();
-<<<<<<< Updated upstream
-=======
         updatePrimaryVariables(simulator, well_state, deferred_logger);
         initPrimaryVariablesEvaluation();
->>>>>>> Stashed changes
         const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
 
         if (this->param_.check_well_operability_iter_)
@@ -1170,7 +1167,7 @@ namespace Opm
                 ws.surface_rates[p] = 0;
             }
             ws.thp = 0;
-            updatePrimaryVariables(simulator, well_state, deferred_logger);
+            //this->updatePrimaryVariables(simulator, well_state, deferred_logger); // can not called becase it is notconst
             return;
         }
 
@@ -1525,7 +1522,7 @@ namespace Opm
                 ws.bhp = controls.bhp_limit;
             }
         }
-        updatePrimaryVariables(simulator, well_state, deferred_logger);
+        //this->updatePrimaryVariables(simulator, well_state, deferred_logger);// can not be called because it is not const
     }
 
     template<typename TypeTag>