diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 6cfb634cb..04b86894b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -95,6 +95,7 @@ list (APPEND TEST_SOURCE_FILES # tests/test_thresholdpressure.cpp tests/test_wellswitchlogger.cpp tests/test_timer.cpp + tests/test_invert.cpp ) list (APPEND TEST_DATA_FILES diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index e2bbfdb78..82797ea17 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -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]; diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index 520e90a16..81037142e 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/opm/autodiff/ISTLSolver.hpp @@ -43,6 +43,138 @@ namespace Dune { +namespace FMatrixHelp { +//! invert 4x4 Matrix without changing the original matrix +template +static inline K invertMatrix (const FieldMatrix &matrix, FieldMatrix &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 &matrix) FMatrixHelp::invertMatrix(A, matrix ); } +//! invert matrix by calling FMatrixHelp::invert +template +static inline void invertMatrix (FieldMatrix &matrix) +{ + FieldMatrix A ( matrix ); + FMatrixHelp::invertMatrix(A, matrix ); +} + //! invert matrix by calling matrix.invert template static inline void invertMatrix (FieldMatrix &matrix) diff --git a/tests/test_invert.cpp b/tests/test_invert.cpp new file mode 100644 index 000000000..b6690663f --- /dev/null +++ b/tests/test_invert.cpp @@ -0,0 +1,77 @@ +/* + Copyright 2017 IRIS AS + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE InvertSpecializationTest +#include +#include + +void checkIdentity(Dune::FieldMatrix M) { + double diag = 0.0; + double offDiag = 0.0; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + if (i == j) + diag += M[i][j]; + else + offDiag += M[i][j]; + } + } + BOOST_CHECK_CLOSE(4, diag, 1e-14); + BOOST_CHECK_SMALL(offDiag, 1e-14); +} + +BOOST_AUTO_TEST_CASE(Invert4x4) +{ + typedef Dune::FieldMatrix BaseType; + BaseType matrix; + BaseType eye; + BaseType inverse; + + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + matrix[i][j] = i + 4*j + 1; + } + } + BaseType matrix_sing (matrix); + // make matrix non-singular + matrix[3][0] = 5; + matrix[0][3] = 14; + + double det = Dune::FMatrixHelp::invertMatrix(matrix, inverse); + BOOST_CHECK_CLOSE(4, det, 1e-14); + + // check matrix * inverse close to identiy + checkIdentity(matrix.rightmultiply(inverse)); + + // check return identity matrix if singular matrix + inverse = 0.0; + double det2 = Dune::FMatrixHelp::invertMatrix(matrix_sing, inverse); + BOOST_CHECK_CLOSE(1.0, det2, 1e-14); + checkIdentity(inverse); +} + + + +