OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-11 00:25:33 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #770 from hnil/matrix_direct_access

Browse Source 
Matrix direct access
This commit is contained in:
Atgeirr Flø Rasmussen 2022-12-19 13:19:20 +01:00 committed by GitHub
commit 29ff1c75f7
2 changed files with 25 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
opm/models/discretization/common/tpfalinearizer.hh
View File

@ -348,7 +348,7 @@ private:
if (materialLawManager->hasDirectionalRelperms()) { if (materialLawManager->hasDirectionalRelperms()) {
dirId = scvf.faceDirFromDirId(); dirId = scvf.faceDirFromDirId();
} }
loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId}; loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId, nullptr};
} }
} }
neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end()); neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
@ -385,9 +385,16 @@ private:
// allocate raw matrix // allocate raw matrix
jacobian_.reset(new SparseMatrixAdapter(simulator_())); jacobian_.reset(new SparseMatrixAdapter(simulator_()));
diagMatAddress_.resize(numCells);
// create matrix structure based on sparsity pattern // create matrix structure based on sparsity pattern
jacobian_->reserve(sparsityPattern); jacobian_->reserve(sparsityPattern);
for (unsigned globI = 0; globI < numCells; globI++) {
const auto& nbInfos = neighborInfo_[globI];
diagMatAddress_[globI] = jacobian_->blockAddress(globI, globI);
for (auto& nbInfo : nbInfos) {
nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI);
}
}
} }
// reset the global linear system of equations. // reset the global linear system of equations.
@ -454,9 +461,11 @@ private:
adres *= nbInfo.faceArea; adres *= nbInfo.faceArea;
setResAndJacobi(res, bMat, adres); setResAndJacobi(res, bMat, adres);
residual_[globI] += res; residual_[globI] += res;
jacobian_->addToBlock(globI, globI, bMat); //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
*diagMatAddress_[globI] += bMat;
bMat *= -1.0; bMat *= -1.0;
jacobian_->addToBlock(globJ, globI, bMat); //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
*nbInfo.matBlockAddress += bMat;
++loc; ++loc;
} }
@ -477,7 +486,8 @@ private:
bMat *= storefac; bMat *= storefac;
// residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac; // residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac;
residual_[globI] += res; residual_[globI] += res;
jacobian_->addToBlock(globI, globI, bMat); //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
*diagMatAddress_[globI] += bMat;
// wells sources for now (should be moved out) // wells sources for now (should be moved out)
if (well_local) { if (well_local) {
res = 0.0; res = 0.0;
@ -487,7 +497,8 @@ private:
adres *= -volume; adres *= -volume;
setResAndJacobi(res, bMat, adres); setResAndJacobi(res, bMat, adres);
residual_[globI] += res; residual_[globI] += res;
jacobian_->addToBlock(globI, globI, bMat); //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
*diagMatAddress_[globI] += bMat;
} }
} // end of loop for cell globI. } // end of loop for cell globI.
@ -505,7 +516,8 @@ private:
adres *= bdyInfo.bcdata.faceArea; adres *= bdyInfo.bcdata.faceArea;
setResAndJacobi(res, bMat, adres); setResAndJacobi(res, bMat, adres);
residual_[globI] += res; residual_[globI] += res;
jacobian_->addToBlock(globI, globI, bMat); ////SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
*diagMatAddress_[globI] += bMat;
} }
} }
@ -547,8 +559,10 @@ private:
double trans; double trans;
double faceArea; double faceArea;
FaceDir::DirEnum faceDirection; FaceDir::DirEnum faceDirection;
MatrixBlock* matBlockAddress;
}; };
SparseTable<NeighborInfo> neighborInfo_; SparseTable<NeighborInfo> neighborInfo_;
std::vector<MatrixBlock*> diagMatAddress_;
using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState; using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
struct BoundaryConditionData struct BoundaryConditionData

4
opm/simulators/linalg/istlsparsematrixadapter.hh
View File

@ -155,6 +155,10 @@ public:
void block(const size_t rowIdx, const size_t colIdx, MatrixBlock& value) const void block(const size_t rowIdx, const size_t colIdx, MatrixBlock& value) const
{ value = (*istlMatrix_)[rowIdx][colIdx]; } { value = (*istlMatrix_)[rowIdx][colIdx]; }
MatrixBlockType* blockAddress(const size_t rowIdx, const size_t colIdx) const
{ return &(*istlMatrix_)[rowIdx][colIdx]; }
/*! /*!
* \brief Set matrix block to given block. * \brief Set matrix block to given block.
*/ */