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

Merge pull request #1257 from totto82/add_invert4x4

Browse Source 
Add invert4x4 specialization for invertMatrix()
This commit is contained in:
Atgeirr Flø Rasmussen
2017-09-04 09:08:32 +02:00
committed by GitHub
parent a162d883bf 210f1edf25
commit 77f3bdaeda
4 changed files with 219 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -702,6 +702,7 @@ namespace Opm {
if (has_solvent_) {
double& ss = priVars[Indices::solventSaturationIdx];
ss -= satScaleFactor * dss;
ss = std::min(std::max(ss, 0.0),1.0);
}
if (has_polymer_) {
double& c = priVars[Indices::polymerConcentrationIdx];

140
opm/autodiff/ISTLSolver.hpp
View File

@@ -43,6 +43,138 @@
namespace Dune
{
namespace FMatrixHelp {
//! invert 4x4 Matrix without changing the original matrix
template <typename K>
static inline K invertMatrix (const FieldMatrix<K,4,4> &matrix, FieldMatrix<K,4,4> &inverse)
{
inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
matrix[1][1] * matrix[2][3] * matrix[3][2] -
matrix[2][1] * matrix[1][2] * matrix[3][3] +
matrix[2][1] * matrix[1][3] * matrix[3][2] +
matrix[3][1] * matrix[1][2] * matrix[2][3] -
matrix[3][1] * matrix[1][3] * matrix[2][2];
inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
matrix[1][0] * matrix[2][3] * matrix[3][2] +
matrix[2][0] * matrix[1][2] * matrix[3][3] -
matrix[2][0] * matrix[1][3] * matrix[3][2] -
matrix[3][0] * matrix[1][2] * matrix[2][3] +
matrix[3][0] * matrix[1][3] * matrix[2][2];
inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
matrix[1][0] * matrix[2][3] * matrix[3][1] -
matrix[2][0] * matrix[1][1] * matrix[3][3] +
matrix[2][0] * matrix[1][3] * matrix[3][1] +
matrix[3][0] * matrix[1][1] * matrix[2][3] -
matrix[3][0] * matrix[1][3] * matrix[2][1];
inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
matrix[1][0] * matrix[2][2] * matrix[3][1] +
matrix[2][0] * matrix[1][1] * matrix[3][2] -
matrix[2][0] * matrix[1][2] * matrix[3][1] -
matrix[3][0] * matrix[1][1] * matrix[2][2] +
matrix[3][0] * matrix[1][2] * matrix[2][1];
inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
matrix[0][1] * matrix[2][3] * matrix[3][2] +
matrix[2][1] * matrix[0][2] * matrix[3][3] -
matrix[2][1] * matrix[0][3] * matrix[3][2] -
matrix[3][1] * matrix[0][2] * matrix[2][3] +
matrix[3][1] * matrix[0][3] * matrix[2][2];
inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
matrix[0][0] * matrix[2][3] * matrix[3][2] -
matrix[2][0] * matrix[0][2] * matrix[3][3] +
matrix[2][0] * matrix[0][3] * matrix[3][2] +
matrix[3][0] * matrix[0][2] * matrix[2][3] -
matrix[3][0] * matrix[0][3] * matrix[2][2];
inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
matrix[0][0] * matrix[2][3] * matrix[3][1] +
matrix[2][0] * matrix[0][1] * matrix[3][3] -
matrix[2][0] * matrix[0][3] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[2][3] +
matrix[3][0] * matrix[0][3] * matrix[2][1];
inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
matrix[0][0] * matrix[2][2] * matrix[3][1] -
matrix[2][0] * matrix[0][1] * matrix[3][2] +
matrix[2][0] * matrix[0][2] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[2][2] -
matrix[3][0] * matrix[0][2] * matrix[2][1];
inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
matrix[0][1] * matrix[1][3] * matrix[3][2] -
matrix[1][1] * matrix[0][2] * matrix[3][3] +
matrix[1][1] * matrix[0][3] * matrix[3][2] +
matrix[3][1] * matrix[0][2] * matrix[1][3] -
matrix[3][1] * matrix[0][3] * matrix[1][2];
inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
matrix[0][0] * matrix[1][3] * matrix[3][2] +
matrix[1][0] * matrix[0][2] * matrix[3][3] -
matrix[1][0] * matrix[0][3] * matrix[3][2] -
matrix[3][0] * matrix[0][2] * matrix[1][3] +
matrix[3][0] * matrix[0][3] * matrix[1][2];
inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
matrix[0][0] * matrix[1][3] * matrix[3][1] -
matrix[1][0] * matrix[0][1] * matrix[3][3] +
matrix[1][0] * matrix[0][3] * matrix[3][1] +
matrix[3][0] * matrix[0][1] * matrix[1][3] -
matrix[3][0] * matrix[0][3] * matrix[1][1];
inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
matrix[0][0] * matrix[1][2] * matrix[3][1] +
matrix[1][0] * matrix[0][1] * matrix[3][2] -
matrix[1][0] * matrix[0][2] * matrix[3][1] -
matrix[3][0] * matrix[0][1] * matrix[1][2] +
matrix[3][0] * matrix[0][2] * matrix[1][1];
inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
matrix[0][1] * matrix[1][3] * matrix[2][2] +
matrix[1][1] * matrix[0][2] * matrix[2][3] -
matrix[1][1] * matrix[0][3] * matrix[2][2] -
matrix[2][1] * matrix[0][2] * matrix[1][3] +
matrix[2][1] * matrix[0][3] * matrix[1][2];
inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
matrix[0][0] * matrix[1][3] * matrix[2][2] -
matrix[1][0] * matrix[0][2] * matrix[2][3] +
matrix[1][0] * matrix[0][3] * matrix[2][2] +
matrix[2][0] * matrix[0][2] * matrix[1][3] -
matrix[2][0] * matrix[0][3] * matrix[1][2];
inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
matrix[0][0] * matrix[1][3] * matrix[2][1] +
matrix[1][0] * matrix[0][1] * matrix[2][3] -
matrix[1][0] * matrix[0][3] * matrix[2][1] -
matrix[2][0] * matrix[0][1] * matrix[1][3] +
matrix[2][0] * matrix[0][3] * matrix[1][1];
inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
matrix[0][0] * matrix[1][2] * matrix[2][1] -
matrix[1][0] * matrix[0][1] * matrix[2][2] +
matrix[1][0] * matrix[0][2] * matrix[2][1] +
matrix[2][0] * matrix[0][1] * matrix[1][2] -
matrix[2][0] * matrix[0][2] * matrix[1][1];
K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
// return identity for singular or nearly singular matrices.
if (std::abs(det) < 1e-40) {
for (int i = 0; i < 4; ++i){
inverse[i][i] = 1.0;
}
return 1.0;
}
K inv_det = 1.0 / det;
inverse *= inv_det;
return det;
}
} // end FMatrixHelp
namespace ISTLUtility {
@@ -70,6 +202,14 @@ static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,4,4> &matrix)
{
FieldMatrix<K,4,4> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling matrix.invert
template <typename K, int n>
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)