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

Improved initial guess of L

Browse Source
This commit is contained in:
tqiu
2021-01-28 13:46:01 +01:00
parent 837a83fc16
commit d222dffcfd
1 changed files with 22 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
opm/simulators/linalg/bda/BILU0.cpp
View File

@@ -169,15 +169,32 @@ namespace bda
// Ut stores diagonal for now // Ut stores diagonal for now
// Ut is actually BSC format // Ut is actually BSC format
std::unique_ptr<BlockedMatrix<bs> > Ut = std::make_unique<BlockedMatrix<bs> >(Umat->Nb, Umat->nnzbs + Umat->Nb); std::unique_ptr<BlockedMatrix<bs> > Ut = std::make_unique<BlockedMatrix<bs> >(Nb, (nnzbs + Nb) / 2);
Lmat->rowPointers[0] = 0; Lmat->rowPointers[0] = 0;
for (int i = 0; i < Nb+1; i++) { for (int i = 0; i < Nb+1; i++) {
Ut->rowPointers[i] = 0; Ut->rowPointers[i] = 0;
} }
Opm::Detail::Inverter<bs> inverter;
// store inverted diagonal
for (int i = 0; i < Nb; i++) {
int iRowStart = LUmat->rowPointers[i];
int iRowEnd = LUmat->rowPointers[i + 1];
// for every block in this row
for (int ij = iRowStart; ij < iRowEnd; ij++) {
int j = LUmat->colIndices[ij];
if (i == j) {
inverter(LUmat->nnzValues + ij * bs * bs, invDiagVals + i * bs * bs);
}
}
}
// initialize initial guess for L: L_A * D
// L_A is strictly lower triangular part of A
// D is inv(diag(A))
int num_blocks_L = 0; int num_blocks_L = 0;
// for every row
for (int i = 0; i < Nb; i++) { for (int i = 0; i < Nb; i++) {
int iRowStart = LUmat->rowPointers[i]; int iRowStart = LUmat->rowPointers[i];
int iRowEnd = LUmat->rowPointers[i + 1]; int iRowEnd = LUmat->rowPointers[i + 1];
@@ -188,7 +205,8 @@ namespace bda
Ut->rowPointers[j+1]++; // actually colPointers, now simply indicates how many blocks this col holds Ut->rowPointers[j+1]++; // actually colPointers, now simply indicates how many blocks this col holds
} else { } else {
Lmat->colIndices[num_blocks_L] = j; Lmat->colIndices[num_blocks_L] = j;
std::copy(LUmat->nnzValues + ij * bs * bs, LUmat->nnzValues + (ij+1) * bs * bs, Lmat->nnzValues + num_blocks_L * bs * bs); // multiply block of L with corresponding diag block
blockMult<bs>(LUmat->nnzValues + ij * bs * bs, invDiagVals + i * bs * bs, Lmat->nnzValues + num_blocks_L * bs * bs);
num_blocks_L++; num_blocks_L++;
} }
} }
@@ -199,7 +217,7 @@ namespace bda
// prefix sum to sum rowsizes into colpointers // prefix sum to sum rowsizes into colpointers
std::partial_sum(Ut->rowPointers, Ut->rowPointers+Nb+1, Ut->rowPointers); std::partial_sum(Ut->rowPointers, Ut->rowPointers+Nb+1, Ut->rowPointers);
// for every row // initialize initial guess for U
for (int i = 0; i < Nb; i++) { for (int i = 0; i < Nb; i++) {
int iRowStart = LUmat->rowPointers[i]; int iRowStart = LUmat->rowPointers[i];
int iRowEnd = LUmat->rowPointers[i + 1]; int iRowEnd = LUmat->rowPointers[i + 1];
@@ -222,7 +240,6 @@ namespace bda
} }
Ut->rowPointers[0] = 0; Ut->rowPointers[0] = 0;
Opm::Detail::Inverter<bs> inverter;
// Utmp is needed for CPU and GPU decomposition, because U is transposed, and reversed after decomposition // Utmp is needed for CPU and GPU decomposition, because U is transposed, and reversed after decomposition
// U will be reversed because it is used with backwards substitution, the last row is used first // U will be reversed because it is used with backwards substitution, the last row is used first