mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
cleaning code
This commit is contained in:
committed by
Atgeirr Flø Rasmussen
parent
5a917dd11e
commit
6c407506a9
@@ -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<Dune::FieldMatrix<double, 1, 1>>;
|
||||
|
||||
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<std::vector<int>> getMaxWellConnections() const
|
||||
{
|
||||
std::vector<std::vector<int>> 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<int> 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<int> 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<std::vector<int>> 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<std::vector<int>> 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<WellInterfacePtr>& localNonshutWells()
|
||||
{
|
||||
return well_container_;
|
||||
}
|
||||
const std::vector<WellInterfacePtr>& localNonshutWells() const;
|
||||
|
||||
int numLocalNonshutWells() const
|
||||
{
|
||||
return well_container_.size();
|
||||
}
|
||||
|
||||
protected:
|
||||
Simulator& ebosSimulator_;
|
||||
|
||||
@@ -1173,10 +1173,140 @@ namespace Opm {
|
||||
Ax.axpy( alpha, scaleAddRes_ );
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilWellModel<TypeTag>::
|
||||
addWellContributions(SparseMatrixAdapter& jacobian) const
|
||||
{
|
||||
for ( const auto& well: well_container_ ) {
|
||||
well->addWellContributions(jacobian);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilWellModel<TypeTag>::
|
||||
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<typename TypeTag>
|
||||
int
|
||||
BlackoilWellModel<TypeTag>::
|
||||
numLocalWellsEnd() const
|
||||
{
|
||||
auto w = schedule().getWellsatEnd();
|
||||
w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end());
|
||||
return w.size();
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
std::vector<std::vector<int>>
|
||||
BlackoilWellModel<TypeTag>::
|
||||
getMaxWellConnections() const
|
||||
{
|
||||
std::vector<std::vector<int>> 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<int> 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<int> 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<typename TypeTag>
|
||||
void
|
||||
BlackoilWellModel<TypeTag>::
|
||||
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<std::vector<int>> 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<typename TypeTag>
|
||||
const std::vector<WellInterfacePtr>&
|
||||
BlackoilWellModel<TypeTag>::
|
||||
localNonshutWells() const;
|
||||
{
|
||||
return well_container_;
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilWellModel<TypeTag>::
|
||||
numLocalNonshutWells() const
|
||||
{
|
||||
return well_container_.size();
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilWellModel<TypeTag>::
|
||||
|
||||
@@ -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);
|
||||
// }
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -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,
|
||||
|
||||
@@ -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 <typename TypeTag>
|
||||
// void
|
||||
// StandardWell<TypeTag>::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 <typename TypeTag>
|
||||
void
|
||||
StandardWell<TypeTag>::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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -232,7 +232,8 @@ namespace Opm {
|
||||
return cube;
|
||||
}
|
||||
|
||||
|
||||
// explicite transpose of dense matrix due to compilation problems
|
||||
// used for caclulating quasiimpes well weights
|
||||
template <class DenseMatrix>
|
||||
DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
|
||||
{
|
||||
|
||||
@@ -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;
|
||||
|
||||
|
||||
@@ -534,7 +534,33 @@ namespace Opm
|
||||
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
|
||||
}
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
WellInterface<TypeTag>::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<typename TypeTag>
|
||||
void
|
||||
|
||||
Reference in New Issue
Block a user