From 888565c335124d59f1034a0983b5d7a1f6f5fb0c Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Dec 2022 11:59:47 +0100 Subject: [PATCH 1/5] start on direct access to matrices --- .../discretization/common/tpfalinearizer.hh | 27 ++++++++++++++----- .../linalg/istlsparsematrixadapter.hh | 4 +++ 2 files changed, 24 insertions(+), 7 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index a5f6161cd..0080cef82 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -347,7 +347,7 @@ private: if (materialLawManager->hasDirectionalRelperms()) { 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()); @@ -384,9 +384,15 @@ private: // allocate raw matrix jacobian_.reset(new SparseMatrixAdapter(simulator_())); - + diagElemAddress_.resize(numCells); // create matrix structure based on sparsity pattern jacobian_->reserve(sparsityPattern); + for (unsigned globI = 0; globI < numCells; globI++) { + diagElements_[globI] = jacobian_->blockAddress(globI,globI,bMat); + for (auto& nbInfo : nbInfos) { + nbInfo.matAdress = jacobian_->blockAddress(globI,nbInfo.neighbor); + } + } } // reset the global linear system of equations. @@ -453,9 +459,11 @@ private: adres *= nbInfo.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - jacobian_->addToBlock(globI, globI, bMat); + //jacobian_->addToBlock(globI, globI, bMat); + *diagElemAddress_[globI] += bMat; bMat *= -1.0; - jacobian_->addToBlock(globJ, globI, bMat); + //jacobian_->addToBlock(globJ, globI, bMat); + *nbInfo.matAddress += bMat; ++loc; } @@ -476,7 +484,8 @@ private: bMat *= storefac; // residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac; residual_[globI] += res; - jacobian_->addToBlock(globI, globI, bMat); + //jacobian_->addToBlock(globI, globI, bMat); + *diagElemAddress_[globI] += bMat; // wells sources for now (should be moved out) if (well_local) { res = 0.0; @@ -486,7 +495,8 @@ private: adres *= -volume; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - jacobian_->addToBlock(globI, globI, bMat); + //jacobian_->addToBlock(globI, globI, bMat); + *diagElemAddress_[globI] += bMat; } } // end of loop for cell globI. @@ -504,7 +514,8 @@ private: adres *= bdyInfo.bcdata.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - jacobian_->addToBlock(globI, globI, bMat); + //jacobian_->addToBlock(globI, globI, bMat); + *diagElemAddress_[globI] += bMat; } } @@ -546,8 +557,10 @@ private: double trans; double faceArea; FaceDir::DirEnum faceDirection; + MatrixBlock* matAddres; }; SparseTable neighborInfo_; + std::vector diagElemAddress_; using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState; struct BoundaryConditionData diff --git a/opm/simulators/linalg/istlsparsematrixadapter.hh b/opm/simulators/linalg/istlsparsematrixadapter.hh index 56582f0ca..a2b8a44d0 100644 --- a/opm/simulators/linalg/istlsparsematrixadapter.hh +++ b/opm/simulators/linalg/istlsparsematrixadapter.hh @@ -155,6 +155,10 @@ public: void block(const size_t rowIdx, const size_t colIdx, MatrixBlock& value) const { 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. */ From d465ebdb7843ffccf0fa71f35147b0eff421df74 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Dec 2022 20:01:35 +0100 Subject: [PATCH 2/5] fixed direct acces to matrix in lineariser --- .../discretization/common/tpfalinearizer.hh | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 0080cef82..d3b0e995d 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -384,13 +384,14 @@ private: // allocate raw matrix jacobian_.reset(new SparseMatrixAdapter(simulator_())); - diagElemAddress_.resize(numCells); + diagMatAddress_.resize(numCells); // create matrix structure based on sparsity pattern jacobian_->reserve(sparsityPattern); for (unsigned globI = 0; globI < numCells; globI++) { - diagElements_[globI] = jacobian_->blockAddress(globI,globI,bMat); + const auto& nbInfos = neighborInfo_[globI]; + diagMatAddress_[globI] = jacobian_->blockAddress(globI,globI); for (auto& nbInfo : nbInfos) { - nbInfo.matAdress = jacobian_->blockAddress(globI,nbInfo.neighbor); + nbInfo.matAddress = jacobian_->blockAddress(nbInfo.neighbor,globI); } } } @@ -435,23 +436,23 @@ private: MatrixBlock bMat(0.0); ADVectorBlock adres(0.0); const IntensiveQuantities* intQuantsInP = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0); - if (intQuantsInP == nullptr) { - throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI)); - } + // if (intQuantsInP == nullptr) { + // throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI)); + // } const IntensiveQuantities& intQuantsIn = *intQuantsInP; // Flux term. short loc = 0; for (const auto& nbInfo : nbInfos) { unsigned globJ = nbInfo.neighbor; - assert(globJ != globI); + //assert(globJ != globI); res = 0.0; bMat = 0.0; adres = 0.0; const IntensiveQuantities* intQuantsExP = model_().cachedIntensiveQuantities(globJ, /*timeIdx*/ 0); - if (intQuantsExP == nullptr) { - throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globJ) + " when assembling fluxes for cell " + std::to_string(globI)); - } + // if (intQuantsExP == nullptr) { + // throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globJ) + " when assembling fluxes for cell " + std::to_string(globI)); + // } const IntensiveQuantities& intQuantsEx = *intQuantsExP; LocalResidual::computeFlux( adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, @@ -460,7 +461,7 @@ private: setResAndJacobi(res, bMat, adres); residual_[globI] += res; //jacobian_->addToBlock(globI, globI, bMat); - *diagElemAddress_[globI] += bMat; + *diagMatAddress_[globI] += bMat; bMat *= -1.0; //jacobian_->addToBlock(globJ, globI, bMat); *nbInfo.matAddress += bMat; @@ -485,7 +486,7 @@ private: // residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac; residual_[globI] += res; //jacobian_->addToBlock(globI, globI, bMat); - *diagElemAddress_[globI] += bMat; + *diagMatAddress_[globI] += bMat; // wells sources for now (should be moved out) if (well_local) { res = 0.0; @@ -496,7 +497,7 @@ private: setResAndJacobi(res, bMat, adres); residual_[globI] += res; //jacobian_->addToBlock(globI, globI, bMat); - *diagElemAddress_[globI] += bMat; + *diagMatAddress_[globI] += bMat; } } // end of loop for cell globI. @@ -515,7 +516,7 @@ private: setResAndJacobi(res, bMat, adres); residual_[globI] += res; //jacobian_->addToBlock(globI, globI, bMat); - *diagElemAddress_[globI] += bMat; + *diagMatAddress_[globI] += bMat; } } @@ -557,10 +558,10 @@ private: double trans; double faceArea; FaceDir::DirEnum faceDirection; - MatrixBlock* matAddres; + MatrixBlock* matAddress; }; SparseTable neighborInfo_; - std::vector diagElemAddress_; + std::vector diagMatAddress_; using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState; struct BoundaryConditionData From 83e7bbb140b20a05ad08310626d9d036b0234eb5 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 19 Dec 2022 11:32:49 +0100 Subject: [PATCH 3/5] addressed review comments --- .../discretization/common/tpfalinearizer.hh | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index d3b0e995d..55f7160c0 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -389,9 +389,9 @@ private: jacobian_->reserve(sparsityPattern); for (unsigned globI = 0; globI < numCells; globI++) { const auto& nbInfos = neighborInfo_[globI]; - diagMatAddress_[globI] = jacobian_->blockAddress(globI,globI); + diagMatAddress_[globI] = jacobian_->blockAddress(globI, globI); for (auto& nbInfo : nbInfos) { - nbInfo.matAddress = jacobian_->blockAddress(nbInfo.neighbor,globI); + nbInfo.matBlockAddress = jacobian_->blockAddress(nbInfo.neighbor, globI); } } } @@ -450,9 +450,9 @@ private: bMat = 0.0; adres = 0.0; const IntensiveQuantities* intQuantsExP = model_().cachedIntensiveQuantities(globJ, /*timeIdx*/ 0); - // if (intQuantsExP == nullptr) { - // throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globJ) + " when assembling fluxes for cell " + std::to_string(globI)); - // } + if (intQuantsExP == nullptr) { + throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globJ) + " when assembling fluxes for cell " + std::to_string(globI)); + } const IntensiveQuantities& intQuantsEx = *intQuantsExP; LocalResidual::computeFlux( adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, @@ -460,11 +460,11 @@ private: adres *= nbInfo.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - //jacobian_->addToBlock(globI, globI, bMat); + //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); *diagMatAddress_[globI] += bMat; bMat *= -1.0; - //jacobian_->addToBlock(globJ, globI, bMat); - *nbInfo.matAddress += bMat; + //SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat); + *nbInfo.matBlockAddress += bMat; ++loc; } @@ -485,7 +485,7 @@ private: bMat *= storefac; // residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac; 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) if (well_local) { @@ -496,7 +496,7 @@ private: adres *= -volume; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - //jacobian_->addToBlock(globI, globI, bMat); + //SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); *diagMatAddress_[globI] += bMat; } } // end of loop for cell globI. @@ -558,7 +558,7 @@ private: double trans; double faceArea; FaceDir::DirEnum faceDirection; - MatrixBlock* matAddress; + MatrixBlock* matBlockAddress; }; SparseTable neighborInfo_; std::vector diagMatAddress_; From e5c76f6eb20f69be49808fb3885ec1c1fd4155de Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 19 Dec 2022 11:34:06 +0100 Subject: [PATCH 4/5] fixed rest of review comments --- opm/models/discretization/common/tpfalinearizer.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 55f7160c0..4063192fc 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -436,9 +436,9 @@ private: MatrixBlock bMat(0.0); ADVectorBlock adres(0.0); const IntensiveQuantities* intQuantsInP = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0); - // if (intQuantsInP == nullptr) { - // throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI)); - // } + if (intQuantsInP == nullptr) { + throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI)); + } const IntensiveQuantities& intQuantsIn = *intQuantsInP; // Flux term. From da1c1747706b269de782a07b099b0c815f0a1357 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 19 Dec 2022 11:38:16 +0100 Subject: [PATCH 5/5] more review comments addressed --- opm/models/discretization/common/tpfalinearizer.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 4063192fc..bef855685 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -445,7 +445,7 @@ private: short loc = 0; for (const auto& nbInfo : nbInfos) { unsigned globJ = nbInfo.neighbor; - //assert(globJ != globI); + assert(globJ != globI); res = 0.0; bMat = 0.0; adres = 0.0; @@ -515,7 +515,7 @@ private: adres *= bdyInfo.bcdata.faceArea; setResAndJacobi(res, bMat, adres); residual_[globI] += res; - //jacobian_->addToBlock(globI, globI, bMat); + ////SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat); *diagMatAddress_[globI] += bMat; } }