Avoid copying of matrices inside StandardWellsDense.

This commit is contained in:
Robert Kloefkorn 2016-09-13 20:40:30 +02:00
parent cd749b3452
commit a9663447c2

View File

@ -1,6 +1,7 @@
/* /*
Copyright 2016 SINTEF ICT, Applied Mathematics. Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA. Copyright 2016 Statoil ASA.
Copyright 2016 IRIS AS
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -95,6 +96,8 @@ namespace Opm {
const std::vector<double>& pv) const std::vector<double>& pv)
: wells_active_(wells_arg!=nullptr) : wells_active_(wells_arg!=nullptr)
, wells_(wells_arg) , wells_(wells_arg)
, param_(param)
, terminal_output_(terminal_output)
, fluid_(nullptr) , fluid_(nullptr)
, active_(nullptr) , active_(nullptr)
, vfp_properties_(nullptr) , vfp_properties_(nullptr)
@ -102,9 +105,10 @@ namespace Opm {
, well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells]) , well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells])
, wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases)
, F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases)
, param_(param)
, terminal_output_(terminal_output)
{ {
invDuneD_.setBuildMode( Mat::row_wise );
duneC_.setBuildMode( Mat::row_wise );
duneB_.setBuildMode( Mat::row_wise );
} }
void init(const BlackoilPropsAdInterface* fluid_arg, void init(const BlackoilPropsAdInterface* fluid_arg,
@ -132,18 +136,24 @@ namespace Opm {
const int nw = wells().number_of_wells; const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw]; const int nperf = wells().well_connpos[nw];
int nc = numCells(); const int nc = numCells();
Mat duneD (nw,nw,nw,Mat::row_wise);
Mat duneC (nw,nc,nperf, Mat::row_wise);
Mat duneB (nw,nc,nperf, Mat::row_wise); // B^T
typedef Mat::CreateIterator Iter; // set invDuneD
for(Iter row=duneD.createbegin(); row!=duneD.createend(); ++row){ invDuneD_.setSize( nw, nw, nw );
// Add nonzeros for diagonal
row.insert(row.index()); // set duneC
duneC_.setSize( nw, nc, nperf );
// set duneB
duneB_.setSize( nw, nc, nperf );
for(auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) {
// Add nonzeros for diagonal
row.insert(row.index());
} }
for(Iter row=duneC.createbegin(); row!=duneC.createend(); ++row){
// Add nonzeros for diagonal for(auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) {
// Add nonzeros for diagonal
for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) {
const int cell_idx = wells().well_cells[perf]; const int cell_idx = wells().well_cells[perf];
row.insert(cell_idx); row.insert(cell_idx);
@ -151,18 +161,14 @@ namespace Opm {
} }
// make the B^T matrix // make the B^T matrix
for(Iter row=duneB.createbegin(); row!=duneB.createend(); ++row){ for(auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) { for (int perf = wells().well_connpos[row.index()] ; perf < wells().well_connpos[row.index()+1]; ++perf) {
const int cell_idx = wells().well_cells[perf]; const int cell_idx = wells().well_cells[perf];
row.insert(cell_idx); row.insert(cell_idx);
} }
} }
invDuneD_ = duneD;
duneC_ = duneC;
duneB_ = duneB;
BVector resWell(nw); resWell_.resize( nw );
resWell_ = resWell;
} }
@ -267,7 +273,7 @@ namespace Opm {
} }
// do the local inversion of D. // do the local inversion of D.
localInvert(invDuneD_); localInvert( invDuneD_ );
} }
@ -1377,12 +1383,11 @@ namespace Opm {
std::vector<EvalWell> wellVariables_; std::vector<EvalWell> wellVariables_;
std::vector<double> F0_; std::vector<double> F0_;
//Mat duneA_;
Mat duneB_; Mat duneB_;
Mat duneC_; Mat duneC_;
Mat invDuneD_; Mat invDuneD_;
BVector resWell_; BVector resWell_;
//BVector rhs_;
// protected methods // protected methods
EvalWell getBhp(const int wellIdx) const { EvalWell getBhp(const int wellIdx) const {