diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index a7ed3b5ae..80e735e86 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -264,13 +264,8 @@ namespace Opm { const SimulatorReportSingle& lastReport() const; - void addWellContributions(SparseMatrixAdapter& jacobian) const - { - for ( const auto& well: well_container_ ) { - well->addWellContributions(jacobian); - } - } - + void addWellContributions(SparseMatrixAdapter& jacobian) const; + // called at the beginning of a report step void beginReportStep(const int time_step); @@ -296,118 +291,20 @@ namespace Opm { using PressureMatrix = Dune::BCRSMatrix>; - int numLocalWellsEnd() const - { - auto w = schedule().getWellsatEnd(); - w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end()); - return w.size(); - } + int numLocalWellsEnd() const; - void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const - { - int nw = this->numLocalWellsEnd(); - int rdofs = local_num_cells_; - for(int i=0; i < nw; i++){ - int wdof = rdofs + i; - jacobian[wdof][wdof] = 1.0;// better scaling ? - } + void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const; - for (const auto& well : well_container_) { - well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState()); - } - } - - std::vector> getMaxWellConnections() const - { - std::vector> wells; - // Create cartesian to compressed mapping - const auto& globalCell = grid().globalCell(); - const auto& cartesianSize = grid().logicalCartesianSize(); - - auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2]; - - std::vector cartesianToCompressed(size, -1); - auto begin = globalCell.begin(); - - for ( auto cell = begin, end= globalCell.end(); cell != end; ++cell ) - { - cartesianToCompressed[ *cell ] = cell - begin; - } - - auto schedule_wells = schedule().getWellsatEnd(); - schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end()); - wells.reserve(schedule_wells.size()); - - // initialize the additional cell connections introduced by wells. - for ( const auto& well : schedule_wells ) - { - std::vector compressed_well_perforations; - // All possible completions of the well - const auto& completionSet = well.getConnections(); - compressed_well_perforations.reserve(completionSet.size()); - - for ( size_t c=0; c < completionSet.size(); c++ ) - { - const auto& completion = completionSet.get(c); - int i = completion.getI(); - int j = completion.getJ(); - int k = completion.getK(); - int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k); - int compressed_idx = cartesianToCompressed[cart_grid_idx]; - - if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells. - { - compressed_well_perforations.push_back(compressed_idx); - } - } - - if ( ! compressed_well_perforations.empty() ) - { - std::sort(compressed_well_perforations.begin(), - compressed_well_perforations.end()); - - wells.push_back(compressed_well_perforations); - } - } - return wells; - } + std::vector> getMaxWellConnections() const; - - void addWellPressureEquationsStruct(PressureMatrix& jacobian) const - { - int nw = this->numLocalWellsEnd(); - int rdofs = local_num_cells_; - for(int i=0; i < nw; i++){ - int wdof = rdofs + i; - jacobian.entry(wdof,wdof) = 1.0;// better scaling ? - } - std::vector> wellconnections = getMaxWellConnections(); - for(int i=0; i < nw; i++){ - const auto& perfcells = wellconnections[i]; - for(int perfcell : perfcells){ - int wdof = rdofs + i; - jacobian.entry(wdof,perfcell) = 0.0; - jacobian.entry(perfcell, wdof) = 0.0; - } - } - // for (const auto& well : well_container_) { - // well->addWellPressureEquationsStruct(jacobian); - // } - } - + void addWellPressureEquationsStruct(PressureMatrix& jacobian) const; void initGliftEclWellMap(GLiftEclWells &ecl_well_map); /// \brief Get list of local nonshut wells - const std::vector& localNonshutWells() - { - return well_container_; - } + const std::vector& localNonshutWells() const; int numLocalNonshutWells() const - { - return well_container_.size(); - } protected: Simulator& ebosSimulator_; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 8dddc8e7f..61e4ee641 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1173,10 +1173,140 @@ namespace Opm { Ax.axpy( alpha, scaleAddRes_ ); } + template + void + BlackoilWellModel:: + addWellContributions(SparseMatrixAdapter& jacobian) const + { + for ( const auto& well: well_container_ ) { + well->addWellContributions(jacobian); + } + } + template + void + BlackoilWellModel:: + addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const + { + int nw = this->numLocalWellsEnd(); + int rdofs = local_num_cells_; + for(int i=0; i < nw; i++){ + int wdof = rdofs + i; + jacobian[wdof][wdof] = 1.0;// better scaling ? + } + + for (const auto& well : well_container_) { + well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState()); + } + } + + template + int + BlackoilWellModel:: + numLocalWellsEnd() const + { + auto w = schedule().getWellsatEnd(); + w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end()); + return w.size(); + } + template + std::vector> + BlackoilWellModel:: + getMaxWellConnections() const + { + std::vector> wells; + // Create cartesian to compressed mapping + const auto& globalCell = grid().globalCell(); + const auto& cartesianSize = grid().logicalCartesianSize(); + auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2]; + + std::vector cartesianToCompressed(size, -1); + auto begin = globalCell.begin(); + + for ( auto cell = begin, end= globalCell.end(); cell != end; ++cell ) + { + cartesianToCompressed[ *cell ] = cell - begin; + } + + auto schedule_wells = schedule().getWellsatEnd(); + schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end()); + wells.reserve(schedule_wells.size()); + + // initialize the additional cell connections introduced by wells. + for ( const auto& well : schedule_wells ) + { + std::vector compressed_well_perforations; + // All possible completions of the well + const auto& completionSet = well.getConnections(); + compressed_well_perforations.reserve(completionSet.size()); + + for ( size_t c=0; c < completionSet.size(); c++ ) + { + const auto& completion = completionSet.get(c); + int i = completion.getI(); + int j = completion.getJ(); + int k = completion.getK(); + int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k); + int compressed_idx = cartesianToCompressed[cart_grid_idx]; + + if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells. + { + compressed_well_perforations.push_back(compressed_idx); + } + } + + if ( ! compressed_well_perforations.empty() ) + { + std::sort(compressed_well_perforations.begin(), + compressed_well_perforations.end()); + + wells.push_back(compressed_well_perforations); + } + } + return wells; + } + + template + void + BlackoilWellModel:: + addWellPressureEquationsStruct(PressureMatrix& jacobian) const + { + int nw = this->numLocalWellsEnd(); + int rdofs = local_num_cells_; + for(int i=0; i < nw; i++){ + int wdof = rdofs + i; + jacobian.entry(wdof,wdof) = 1.0;// better scaling ? + } + std::vector> wellconnections = getMaxWellConnections(); + for(int i=0; i < nw; i++){ + const auto& perfcells = wellconnections[i]; + for(int perfcell : perfcells){ + int wdof = rdofs + i; + jacobian.entry(wdof,perfcell) = 0.0; + jacobian.entry(perfcell, wdof) = 0.0; + } + } + } + + template + const std::vector& + BlackoilWellModel:: + localNonshutWells() const; + { + return well_container_; + } + + template + int + BlackoilWellModel:: + numLocalNonshutWells() const + { + return well_container_.size(); + } + template void BlackoilWellModel:: diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index f04354aef..84ff79d22 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -756,96 +756,62 @@ namespace Opm 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. + // Add the pressure contribution to the cpr system for the well + + // Add for coupling from well to reservoir 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))){ + if(not(this->isPressureControlled(well_state))){ + 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; + 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; + // make cpr weights for well by pure avarage of reservoir weights of the perforations + if(not(this->isPressureControlled(well_state))){ + auto well_weight = weights[0]*0.0; + int num_perfs = 0; + for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) { + 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]; + well_weight += bw; + num_perfs += 1; + } } - //segment_weights[rowB] /= num_perfs; - } - well_weight /= num_perfs; - assert(num_perfs>0); - - double diag_ell = 0.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))){ + + well_weight /= num_perfs; + assert(num_perfs>0); + + // Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well + double diag_ell = 0.0; + for (size_t rowB = 0; rowB < this->duneB_.N(); ++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; for(int i = 0; i< bw.size(); ++i){ matel += bw[i] *(*colB)[i][pressureVarIndex]; } jacobian[welldof_ind][col_index] += matel; diag_ell -= matel; - } + } } - } - - if(this->isPressureControlled(well_state)){ - jacobian[welldof_ind][welldof_ind] = 1.0; - }else{ assert(diag_ell > 0.0); jacobian[welldof_ind][welldof_ind] = diag_ell; + }else{ + jacobian[welldof_ind][welldof_ind] = 1.0; // maybe we could have used diag_ell if calculated } - // 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); - // } - - } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 8389b9c6e..b1ad82a92 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -182,8 +182,6 @@ namespace Opm virtual void addWellContributions(SparseMatrixAdapter& mat) const override; - // virtual void addWellPressureEquationsStruct(PressureMatrix& mat) const override; - virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x, const int pressureVarIndex, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 36de07a69..20d4a66b5 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -550,7 +550,7 @@ namespace Opm // do the local inversion of D. try { - this->duneD_ = this->invDuneD_; + this->duneD_ = this->invDuneD_; // Not strictly need if not cpr with well contributions is used Dune::ISTLUtility::invertMatrix(this->invDuneD_[0][0]); } catch( ... ) { OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger); @@ -2167,39 +2167,6 @@ namespace Opm } } - - - - - // template - // void - // StandardWell::addWellPressureEquationsStruct(PressureMatrix& jacobian) const - // { - // // sustem is the pressur variant of - // // We need to change matrx A as follows - // // A CT - // // B D - // // 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 - // // BHP contolled - // 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_index = colC.index(); - // double matel = 0; - // 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; - // } - // } - - template void StandardWell::addWellPressureEquations(PressureMatrix& jacobian, @@ -2208,76 +2175,65 @@ namespace Opm const bool use_well_weights, const WellState& well_state) const { - // sustem is the pressur variant of - // We need to change matrx A as follows - // A CT - // B D - // 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 - // BHP contolled + // Add the well contributions in cpr + // use_well_weights is a quasiimpes formulation which is not implemented in multisegment int bhp_var_index = Bhp; int nperf = 0; - auto cell_weights = weights[0]*0.0; + auto cell_weights = weights[0]*0.0;// not need for not(use_well_weights) 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; - if(not(this->isPressureControlled(well_state))){ + // do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also) + if(not(this->isPressureControlled(well_state)) || use_well_weights){ + // make coupling for reservoir to 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) { + const auto row_ind = colC.index(); + const auto& bw = weights[row_ind]; + double matel = 0; assert((*colC).M() == bw.size()); for (size_t i = 0; i < bw.size(); ++i) { matel += (*colC)[bhp_var_index][i] * bw[i]; } + + jacobian[row_ind][welldof_ind] = matel; + cell_weights += bw; + nperf += 1; } - 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]; - BVectorWell rhs(1); - rhs[0].resize(blockSz); - rhs[0][bhp_var_index] = 1.0; - DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0]; - DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block); - // diag_block_transpose.solve(bweights, rhs); - //HACK due to template errors + { if(use_well_weights){ + // calculate weighs and set diagonal element + //NB! use this options without treating pressure controlled separated + //NB! calculate quasiimpes well weights NB do not work well with trueimpes reservoir weights double abs_max = 0; - if(this->isPressureControlled(well_state)){ - // examples run ok without this branch also - 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])); + BVectorWell rhs(1); + rhs[0].resize(blockSz); + rhs[0][bhp_var_index] = 1.0; + DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0]; + DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block); + 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]; } - 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])); } + assert(abs_max>0.0); + for (size_t i = 0; i < blockSz; ++i) { + bweights[0][i] /= abs_max; + } + diagElem = 1.0/abs_max; }else{ + // set diagonal element if(this->isPressureControlled(well_state)){ bweights[0][blockSz-1] = 1.0; - diagElem = 1.0;// better scaling? + diagElem = 1.0;// better scaling could have used the calculation below if weights were calculated }else{ for (size_t i = 0; i < cell_weights.size(); ++i) { bweights[0][i] = cell_weights[i]; @@ -2291,26 +2247,20 @@ namespace Opm } } - //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); }); - - - } // jacobian[welldof_ind][welldof_ind] = diagElem; // 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(); - const auto& bw = bweights[0]; - double matel = 0; - for (size_t i = 0; i < bw.size(); ++i) { - const double w = bw[i]; - matel += (*colB)[i][pressureVarIndex] * bw[i]; - } - jacobian[welldof_ind][col_index] = matel; + if(not(this->isPressureControlled(well_state)) || use_well_weights){ + for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) { + const auto col_index = colB.index(); + const auto& bw = bweights[0]; + double matel = 0; + for (size_t i = 0; i < bw.size(); ++i) { + const double w = bw[i]; + matel += (*colB)[i][pressureVarIndex] * bw[i]; + } + jacobian[welldof_ind][col_index] = matel; } } diff --git a/opm/simulators/wells/WellHelpers.hpp b/opm/simulators/wells/WellHelpers.hpp index 9240ec540..e942476a1 100644 --- a/opm/simulators/wells/WellHelpers.hpp +++ b/opm/simulators/wells/WellHelpers.hpp @@ -232,7 +232,8 @@ namespace Opm { return cube; } - + // explicite transpose of dense matrix due to compilation problems + // used for caclulating quasiimpes well weights template DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M) { diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 18f75ab03..efe9b5978 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -225,47 +225,13 @@ public: // Add well contributions to matrix virtual void addWellContributions(SparseMatrixAdapter&) const = 0; - 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 bool isPressureControlled(const WellState& well_state) const; - - // virtual void addWellPressureEquationsStruct(PressureMatrix&) const - // { - // THROW(std::logic_error, "Not implemented for this welltype "); - // } - virtual void addWellPressureEquations(PressureMatrix& mat, 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; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 093204df7..471bab985 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -534,7 +534,33 @@ namespace Opm assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); } - + template + void + WellInterface::isPressureControlled(const WellState& well_state) const + { + 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; + } template void