small bug fixes.

opm/autodiff/StandardWell.hpp

@@ -57,6 +57,12 @@ namespace Opm
         using Base::has_polymer;
         using Base::has_energy;
 
+        // polymer concentration and temperature are already known by the well, so
+        // polymer and energy conservation do not need to be considered explicitly
+        static const int numPolymerEq = has_polymer ? 1 : 0;
+        static const int numEnergyEq = has_energy ? 1 : 0;
+        static const int numWellEq = numEq + 1 - numPolymerEq - numEnergyEq;
+
         // the positions of the primary variables for StandardWell
         // there are four primary variables, the second and the third ones are F_w and F_g
         // the first one is the weighted total rate (G_t), the second and the third ones are F_w and F_g
@@ -64,17 +70,19 @@ namespace Opm
         // the fraction of the solvent, as an extension of the blackoil model, is behind the BHP
         // correspondingly, we have four well equations for blackoil model, the first three are mass
         // converstation equations, and the last one is the well control equation.
+        // primary variables related to other components, will be before the Bhp and after F_g.
+        // well control equation is always the last well equation, other equations will be before the
+        // well control equation and are conservation equations for components involved.
         // TODO: in the current implementation, we use the well rate as the first primary variables for injectors
         // TODO: not sure we should change it.
-
         static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
         static const int GTotal = 0;
         static const int WFrac = gasoil? -1000: 1;
         static const int GFrac = gasoil? 1: 2;
-        // TODO: it is possible the order of Bhp and SFrac need to switched, due to scalingFactor function
+        static const int SFrac = !has_solvent ? -1000 : 3;
-        // TODO: we will do that when we see the problem.
+        // the index for Bhp in primary variables and also the index of well control equation
-        static const int Bhp = gasoil? 2 : 3;
+        // they both will be the last one in their system.
-        static const int SFrac = !has_solvent ? -1000 : Bhp + 1;
+        static const int Bhp = numWellEq - 1;
 
         using typename Base::Scalar;
         using typename Base::ConvergenceReport;
@@ -85,11 +93,6 @@ namespace Opm
         using Base::Oil;
         using Base::Gas;
 
-        // polymer concentration and temperature are already known by the well, so
-        // polymer and energy conservation do not need to be considered explicitly
-        static const int numPolymerEq = has_polymer ? 1 : 0;
-        static const int numEnergyEq = has_energy ? 1 : 0;
-        static const int numWellEq =numEq - numPolymerEq - numEnergyEq;
         using typename Base::Mat;
         using typename Base::BVector;
         using typename Base::Eval;
@@ -239,7 +242,7 @@ namespace Opm
         // TODO: this function should be moved to the base class.
         // while it faces chanllenges for MSWell later, since the calculation of bhp
         // based on THP is never implemented for MSWell yet.
-        EvalWell getBhp() const;
+        const EvalWell& getBhp() const;
 
         // TODO: it is also possible to be moved to the base class.
         EvalWell getQs(const int comp_idx) const;

opm/autodiff/StandardWell_impl.hpp

@@ -105,10 +105,7 @@ namespace Opm
     void StandardWell<TypeTag>::
     initPrimaryVariablesEvaluation() const
     {
-        // TODO: using num_components_ here is only to make the 2p + dummy phase work
+        for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
-        // TODO: in theory, we should use numWellEq here.
-        // for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
-        for (int eqIdx = 0; eqIdx < num_components_; ++eqIdx) {
             assert( (size_t)eqIdx < primary_variables_.size() );
 
             primary_variables_evaluation_[eqIdx] = 0.0;
@@ -122,7 +119,7 @@ namespace Opm
 
 
     template<typename TypeTag>
-    typename StandardWell<TypeTag>::EvalWell
+    const typename StandardWell<TypeTag>::EvalWell&
     StandardWell<TypeTag>::
     getBhp() const
     {
@@ -148,10 +145,10 @@ namespace Opm
     template<typename TypeTag>
     typename StandardWell<TypeTag>::EvalWell
     StandardWell<TypeTag>::
-    getQs(const int comp_idx) const // TODO: phase or component?
+    getQs(const int comp_idx) const
     {
         // TODO: not sure the best way to handle solvent injection
-        // TODO: we need to come back to hanlde the solvent case here, the following implementation does not
+        // TODO: we need to come back to handle the solvent case here, the following implementation does not
         // TODO: consider solvent injection yet.
 
         // TODO: currently, the GTotal definition is still depends on Injector/Producer.
@@ -161,7 +158,6 @@ namespace Opm
             // TODO: using comp_frac here is dangerous, it should be changed later
             // Most likely, it should be changed to use distr, or at least, we need to update comp_frac_ based on distr
             // while solvent might complicate the situation
-            EvalWell qs = 0.0;
             const auto pu = phaseUsage();
             const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx);
             double comp_frac = 0.0;
@@ -746,7 +742,7 @@ namespace Opm
             const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
             const auto& fs = intQuants.fluidState();
             const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
-            const EvalWell bhp = getBhp();
+            const EvalWell& bhp = getBhp();
 
             // Pressure drawdown (also used to determine direction of flow)
             const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
@@ -1045,7 +1041,6 @@ namespace Opm
         // TODO: we should only maintain one current control either from the well_state or from well_controls struct.
         // Either one can be more favored depending on the final strategy for the initilzation of the well control
         const int current = well_state.currentControls()[index_of_well_];
-        const double target_rate = well_controls_iget_target(wc, current);
         const int nwc = well_controls_get_num(wc);
         // Looping over all controls until we find a THP constraint
         for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {