Changes to make ms wells work with cprw

opm/simulators/linalg/PressureBhpTransferPolicy.hpp

@@ -184,7 +184,8 @@ namespace Opm
         }
         if (prm_.get<bool>("add_wells")) {
             assert(transpose == false); // not implemented
-            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_);
+            bool use_well_weights = prm_.get<bool>("use_well_weights");
+            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
 #ifndef NDEBUG
             std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
             assert(rowCoarse == coarseLevelMatrix_->end());

opm/simulators/linalg/WellOperators.hpp

@@ -50,7 +50,7 @@ class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
 {
 public:
     using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
-    virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const = 0;
+    virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const = 0;
     virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
     virtual int getNumberOfExtraEquations() const = 0;
 };
@@ -90,9 +90,9 @@ public:
     {
         return Dune::SolverCategory::sequential;
     }
-    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const override
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const override
     {
-        wellMod_.addWellPressureEquations(jacobian, weights);
+        wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
     }
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
     {
@@ -179,9 +179,9 @@ public:
 
   virtual const matrix_type& getmat() const override { return A_; }
 
-    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
     {
-        wellOper_.addWellPressureEquations(jacobian, weights);
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
     }
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
     {
@@ -269,9 +269,9 @@ public:
 
     virtual const matrix_type& getmat() const override { return A_; }
 
-    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const
+    void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
     {
-        wellOper_.addWellPressureEquations(jacobian, weights);
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
     }
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
     {

opm/simulators/wells/BlackoilWellModel.hpp

@@ -303,7 +303,7 @@ namespace Opm {
                 return w.size();
             }
                    
-            void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
+            void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const
             {
                 int nw =  this->numLocalWellsEnd();
                 int rdofs = local_num_cells_;
@@ -313,7 +313,7 @@ namespace Opm {
                 }
 
                 for (const auto& well : well_container_) {
-                    well->addWellPressureEquations(jacobian, weights, pressureVarIndex);
+                    well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
                 }
             }
 
@@ -390,9 +390,9 @@ namespace Opm {
                         jacobian.entry(perfcell, wdof) = 0.0;
                     }
                 }
-                for (const auto& well : well_container_) {
+                // for (const auto& well : well_container_) {
-                    well->addWellPressureEquationsStruct(jacobian);
+                //     well->addWellPressureEquationsStruct(jacobian);
-                }
+                // }
             }
 
 

opm/simulators/wells/MultisegmentWell.hpp

@@ -74,7 +74,8 @@ namespace Opm
         using MSWEval::GTotal;
         using MSWEval::SPres;
         using MSWEval::numWellEq;
-
+        using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+        
         MultisegmentWell(const Well& well,
                          const ParallelWellInfo& pw_info,
                          const int time_step,
@@ -137,6 +138,12 @@ namespace Opm
                                              DeferredLogger& deferred_logger) const override;
 
         virtual void  addWellContributions(SparseMatrixAdapter& jacobian) const override;
+        
+        virtual void addWellPressureEquations(PressureMatrix& mat,
+                                              const BVector& x,
+                                              const int pressureVarIndex,
+                                              const bool use_well_weights,
+                                              const WellState& well_state) const override;
 
         virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
                                                             DeferredLogger& deferred_logger) const override;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -747,6 +747,103 @@ namespace Opm
     }
 
 
+   template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    addWellPressureEquations(PressureMatrix& jacobian,
+                             const BVector& weights,
+                             const int pressureVarIndex,
+                             const bool use_well_weights,
+                             const WellState& well_state) const
+    {
+        // We need to change matrix A as follows
+        // A -= C^T D^-1 B
+        // D is a (nseg x nseg) block matrix with (4 x 4) blocks.
+        // B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
+        // They have nonzeros at (i, j) only if this well has a
+        // perforation at cell j connected to segment i.  The code
+        // assumes that no cell is connected to more than one segment,
+        // i.e. the columns of B/C have no more than one nonzero.
+        const auto seg_pressure_var_ind  = this->SPres; 
+        const int welldof_ind = this->duneC_.M() + this->index_of_well_;
+        for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
+            for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
+                const auto row_index = colC.index();
+                const auto& bw = weights[row_index];
+                double matel = 0.0;
+                //if(this->isPressureControlled(well_state)){
+                //    jacobian[row_index][welldof_ind] = 0.0;
+                if(not(this->isPressureControlled(well_state))){
+                    for(int i = 0; i< bw.size(); ++i){
+                        matel += bw[i]*(*colC)[seg_pressure_var_ind][i];
+                    }
+                    jacobian[row_index][welldof_ind] += matel;
+                }
+            }
+        }
+
+        //BVector segment_weights(this->duneB_.N());
+        auto well_weight = weights[0]*0.0;
+        int num_perfs = 0;
+        //segment_weights = 0.0;
+        for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
+            //segment_weights[rowB] = 0.0;
+            //int num_perfs = 0;
+            for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
+                const auto col_index = colB.index();
+                const auto& bw = weights[col_index];
+                //segment_weights[rowB] += bw;
+                well_weight += bw;
+                num_perfs += 1;
+            }
+            //segment_weights[rowB] /= num_perfs;
+        }
+        well_weight /= num_perfs;
+        assert(num_perfs>0);
+
+        
+        for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
+            //const auto& bw = segment_weights[rowB];
+            const auto& bw = well_weight;
+            for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
+                const auto col_index = colB.index();               
+                double matel = 0.0;
+                //if(this->isPressureControlled(well_state)){
+                //    jacobian[welldof_ind][col_index] = 0.0;
+                if(not(this->isPressureControlled(well_state))){
+                    for(int i = 0; i< bw.size(); ++i){
+                        matel += bw[i] *(*colB)[i][pressureVarIndex];
+                    }
+                    jacobian[welldof_ind][col_index] += matel;
+                }
+            }
+        }
+
+        
+        if(this->isPressureControlled(well_state)){
+            jacobian[welldof_ind][welldof_ind] = 1.0;
+        }else{
+            for (size_t rowD = 0; rowD < this->duneD_.N(); ++rowD) {
+                //const auto& bw = segment_weights[rowD];
+                const auto& bw = well_weight;
+                //const auto row_index = rowD.index();
+                for (auto colD = this->duneD_[rowD].begin(), endD = this->duneD_[rowD].end(); colD != endD; ++colD) {
+                    const auto col_index = colD.index();
+                    if(rowD == col_index){//need?
+                        double matel = 0.0;
+                        for(int i = 0; i< bw.size(); ++i){
+                            matel += bw[i]*(*colD)[i][seg_pressure_var_ind];
+                        }
+                        jacobian[welldof_ind][welldof_ind] += matel;
+                    }
+                }
+            }
+            assert(jacobian[welldof_ind][welldof_ind] != 0);
+        }
+        
+        
+    }
+    
 
     template<typename TypeTag>
     template<class Value>

opm/simulators/wells/StandardWell.hpp

@@ -182,9 +182,13 @@ namespace Opm
 
         virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
 
-        virtual void addWellPressureEquationsStruct(PressureMatrix& mat) const override;
+        // virtual void addWellPressureEquationsStruct(PressureMatrix& mat) const override;
 
-        virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x,const int pressureVarIndex) const override;
+        virtual void addWellPressureEquations(PressureMatrix& mat,
+                                              const BVector& x,
+                                              const int pressureVarIndex,
+                                              const bool use_well_weights,
+                                              const WellState& well_state) const override;
 
         // iterate well equations with the specified control until converged
         bool iterateWellEqWithControl(const Simulator& ebosSimulator,

opm/simulators/wells/StandardWellGeneric.hpp

@@ -124,6 +124,7 @@ protected:
     OffDiagMatWell duneC_;
     // diagonal matrix for the well
     DiagMatWell invDuneD_;
+    DiagMatWell duneD_;
 
     // Wrapper for the parallel application of B for distributed wells
     wellhelpers::ParallelStandardWellB<Scalar> parallelB_;

opm/simulators/wells/StandardWell_impl.hpp

@@ -550,6 +550,7 @@ namespace Opm
 
         // do the local inversion of D.
         try {
+            this->duneD_ = this->invDuneD_;
             Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]);
         } catch( ... ) {
             OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
@@ -2170,40 +2171,42 @@ namespace Opm
 
 
 
-    template <typename TypeTag>
+    // template <typename TypeTag>
-    void
+    // void
-    StandardWell<TypeTag>::addWellPressureEquationsStruct(PressureMatrix& jacobian) const
+    // StandardWell<TypeTag>::addWellPressureEquationsStruct(PressureMatrix& jacobian) const
-    {
+    // {
-        // sustem is the pressur variant of
+    //     // sustem is the pressur variant of
-        // We need to change matrx A as follows
+    //     // We need to change matrx A as follows
-        // A CT
+    //     // A CT
-        // B D
+    //     // B D
-        // we need to add the elemenst of CT
+    //     // we need to add the elemenst of CT
-        // then we need to ad the quasiimpes type well equation for B D if the well is not
+    //     // then we need to ad the quasiimpes type well equation for B D if the well is not
-        // BHP contolled
+    //     // BHP contolled
-        const int welldof_ind = this->duneC_.M() + this->index_of_well_;
+    //     const int welldof_ind = this->duneC_.M() + this->index_of_well_;
-        for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
+    //     for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
-            const auto row_index = colC.index();
+    //         const auto row_index = colC.index();
-            double matel = 0;
+    //         double matel = 0;
-            jacobian.entry(row_index, welldof_ind) = matel;
+    //         jacobian.entry(row_index, welldof_ind) = matel;
-        }
+    //     }
-
-        jacobian.entry(welldof_ind, welldof_ind) = 0.0;
-
-        // set the matrix elements for well reservoir coupling
-        for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
-            const auto col_index = colB.index();
-            double matel = 0;
-            jacobian.entry(welldof_ind, col_index) = matel;
-        }
-    }
 
+    //     jacobian.entry(welldof_ind, welldof_ind) = 0.0;
 
+    //     // set the matrix elements for well reservoir coupling
+    //     for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
+    //         const auto col_index = colB.index();
+    //         double matel = 0;
+    //         jacobian.entry(welldof_ind, col_index) = matel;
+    //     }
+    // }
 
 
     template <typename TypeTag>
     void
-    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const int pressureVarIndex) const
+    StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
+                                                    const BVector& weights,
+                                                    const int pressureVarIndex,
+                                                    const bool use_well_weights,
+                                                    const WellState& well_state) const
     {
         // sustem is the pressur variant of
         // We need to change matrx A as follows
@@ -2213,24 +2216,34 @@ namespace Opm
         // then we need to ad the quasiimpes type well equation for B D if the well is not
         // BHP contolled
         int bhp_var_index = Bhp;
+        int nperf = 0;
+        auto cell_weights = weights[0]*0.0;
         assert(this->duneC_.M() == weights.size());
         const int welldof_ind = this->duneC_.M() + this->index_of_well_;
         for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
             const auto row_ind = colC.index();
             const auto& bw = weights[row_ind];
             double matel = 0;
-            assert((*colC).M() == bw.size());
+            if(not(this->isPressureControlled(well_state))){
-            for (size_t i = 0; i < bw.size(); ++i) {
+                assert((*colC).M() == bw.size());
-                matel += (*colC)[bhp_var_index][i] * bw[i];
+                for (size_t i = 0; i < bw.size(); ++i) {
+                    matel += (*colC)[bhp_var_index][i] * bw[i];
+                }
             }
             jacobian[row_ind][welldof_ind] = matel;
+            //if(not(use_well_weights)){
+            cell_weights += bw;
+            nperf += 1;
+                //}
         }
+        cell_weights /= nperf;
         // make quasipes weights for bhp it should be trival
         //using VectorBlockType = BVectorWell;
         //VectorBlockType
         BVectorWell  bweights(1);
         size_t blockSz = this->numWellEq_;
         bweights[0].resize(blockSz);
+        bweights[0] = 0.0;
         double diagElem = 0;
         {
             // const DiagMatrixBlockWellType& invA = invDuneD_[0][0];
@@ -2241,23 +2254,49 @@ namespace Opm
             DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
             // diag_block_transpose.solve(bweights, rhs);
             //HACK due to template errors
-            double abs_max = 0;
+            if(use_well_weights){
-            for (size_t i = 0; i < blockSz; ++i) {
+                double abs_max = 0;
-                bweights[0][i] = 0;
+                if(this->isPressureControlled(well_state)){
-                for (size_t j = 0; j < blockSz; ++j) {
+                    // examples run ok without this branch also
-                    bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
+                    bweights[0][blockSz-1] = 1.0;
+                    diagElem = 1.0;// better scaling 
+                }else{
+                    for (size_t i = 0; i < blockSz; ++i) {
+                        bweights[0][i] = 0;
+                        for (size_t j = 0; j < blockSz; ++j) {
+                            bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
+                        }
+                        abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
+                    }
+                    assert(abs_max>0.0);
+                    for (size_t i = 0; i < blockSz; ++i) {
+                        bweights[0][i] /= abs_max;
+                    }
+                    diagElem = 1.0/abs_max;
                 }
-                abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
+            }else{
+                if(this->isPressureControlled(well_state)){
+                    bweights[0][blockSz-1] = 1.0;
+                    diagElem = 1.0;// better scaling?
+                }else{
+                    for (size_t i = 0; i < cell_weights.size(); ++i) {
+                        bweights[0][i] = cell_weights[i];
+                    }
+                    bweights[0][blockSz-1] = 0.0;
+                    diagElem = 0.0;
+                    const auto& locmat = this->duneD_[0][0]; 
+                    for (size_t i = 0; i < cell_weights.size(); ++i) {
+                        diagElem += locmat[i][bhp_var_index]*cell_weights[i];
+                    }
+                    
+                }                 
             }
             //inv_diag_block_transpose.mv(rhs[0], bweights[0]);
             // NB how to scale to make it most symmetric
             //double abs_max = *std::max_element(
             //    bweights[0].begin(), bweights[0].end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
-            assert(abs_max>0.0);
+            
-            for (size_t i = 0; i < blockSz; ++i) {
+            
-                bweights[0][i] /= abs_max;
-            }
-            diagElem = 1.0/abs_max;
             
         }
         //

opm/simulators/wells/WellInterface.hpp

@@ -225,14 +225,47 @@ public:
     // Add well contributions to matrix
     virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
 
-    virtual void addWellPressureEquationsStruct(PressureMatrix&) const
+    virtual bool isPressureControlled(const WellState& well_state) const
     {
+        //return false;
+        bool thp_controlled_well = false;
+        bool bhp_controlled_well = false;
+        const auto& ws = well_state.well(this->index_of_well_);
+        if (this->isInjector()) {
+            const Well::InjectorCMode& current = ws.injection_cmode;
+            if (current == Well::InjectorCMode::THP) {
+                thp_controlled_well = true;
+            }
+            if (current == Well::InjectorCMode::BHP) {
+                bhp_controlled_well = true;
+            }
+        } else {
+            const Well::ProducerCMode& current = ws.production_cmode;
+            if (current == Well::ProducerCMode::THP) {
+                thp_controlled_well = true;
+            }
+            if (current == Well::ProducerCMode::BHP) {
+                bhp_controlled_well = true;
+            }
+        }
+        bool ispressureControlled =  (bhp_controlled_well || thp_controlled_well);
+        return ispressureControlled;
     }
+    
+    
+    // virtual void addWellPressureEquationsStruct(PressureMatrix&) const
+    // {
+    //     THROW(std::logic_error, "Not implemented for this welltype ");
+    // }
 
-    virtual void addWellPressureEquations(PressureMatrix&,
+    virtual void addWellPressureEquations(PressureMatrix& mat,
-                                          const BVector& x,const int pressureVarIndex) const
+                                          const BVector& x,
-    {
+                                          const int pressureVarIndex,
-    }
+                                          const bool use_well_weights,
+                                          const WellState& well_state) const = 0;
+//    {
+        //THROW(std::logic_error, "Not implemented for this welltype ");
+//    }
 
     void addCellRates(RateVector& rates, int cellIdx) const;