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

Fix: zeros were not replaced correctly in the diagonal

Browse Source 
And minor changes
This commit is contained in:
Tong Dong Qiu 2022-06-07 10:21:38 +02:00
parent 26a9582ece
commit 91a54645fb
2 changed files with 9 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
opm/simulators/linalg/bda/BdaBridge.cpp
View File

@ -128,17 +128,16 @@ BdaBridge<BridgeMatrix, BridgeVector, block_size>::BdaBridge(std::string acceler
template <class BridgeMatrix>
int checkZeroDiagonal(BridgeMatrix& mat) {
static std::vector<typename BridgeMatrix::size_type> diag_indices; // contains offsets of the diagonal nnzs
int replaceZeroDiagonal(BridgeMatrix& mat, std::vector<typename BridgeMatrix::size_type>& diag_indices) {
int numZeros = 0;
const int dim = 3; // might be replaced with mat[0][0].N() or BridgeMatrix::block_type::size()
const int dim = mat[0][0].N(); // might be replaced with BridgeMatrix::block_type::size()
const double zero_replace = 1e-15;
if (diag_indices.empty()) {
int N = mat.N();
diag_indices.reserve(N);
int Nb = mat.N();
diag_indices.reserve(Nb);
for (typename BridgeMatrix::iterator r = mat.begin(); r != mat.end(); ++r) {
auto diag = r->find(r.index()); // diag is an iterator
assert(diag.index() == r.index());
assert(diag.index() == r.index()); // every row must have a diagonal block
for (int rr = 0; rr < dim; ++rr) {
auto& val = (*diag)[rr][rr]; // reference to easily change the value
if (val == 0.0) { // could be replaced by '< 1e-30' or similar
@ -244,7 +243,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system(BridgeMatri
}
Dune::Timer t_zeros;
int numZeros = checkZeroDiagonal(*bridgeMat);
int numZeros = replaceZeroDiagonal(*bridgeMat, diagIndices);
if (verbosity >= 2) {
std::ostringstream out;
out << "Checking zeros took: " << t_zeros.stop() << " s, found " << numZeros << " zeros";
@ -263,7 +262,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system(BridgeMatri
}
Dune::Timer t_zeros2;
int jacNumZeros = checkZeroDiagonal(*jacMat);
int jacNumZeros = replaceZeroDiagonal(*jacMat, jacDiagIndices);
if (verbosity >= 2) {
std::ostringstream out;
out << "Checking zeros for jacMat took: " << t_zeros2.stop() << " s, found " << jacNumZeros << " zeros";

2
opm/simulators/linalg/bda/BdaBridge.hpp
View File

@ -48,6 +48,8 @@ private:
std::shared_ptr<Opm::Accelerator::BlockedMatrix> jacMatrix; // 'stores' preconditioner matrix, actually points to h_rows, h_cols and the received BridgeMatrix for the nonzeroes
std::vector<int> h_rows, h_cols; // store the sparsity pattern of the matrix
std::vector<int> h_jacRows, h_jacCols; // store the sparsity pattern of the jacMatrix
std::vector<typename BridgeMatrix::size_type> diagIndices; // contains offsets of the diagonal blocks wrt start of the row, used for replaceZeroDiagonal()
std::vector<typename BridgeMatrix::size_type> jacDiagIndices; // same but for jacMatrix
public:
/// Construct a BdaBridge