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 #337 from atgeirr/improve-helpers

Browse Source 
Add new vertcatCollapseJacs() helper and use it
This commit is contained in:
Bård Skaflestad
2015-03-24 17:13:47 +01:00
parent f288719e07 059bfaf036
commit 7fa2d53db5
4 changed files with 240 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -43,6 +43,7 @@ list (APPEND MAIN_SOURCE_FILES
# originally generated with the command:
# find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
list (APPEND TEST_SOURCE_FILES
tests/test_autodiffhelpers.cpp
tests/test_block.cpp
tests/test_boprops_ad.cpp
tests/test_rateconverter.cpp

92
opm/autodiff/AutoDiffHelpers.hpp
View File

@@ -434,6 +434,98 @@ vertcat(const AutoDiffBlock<double>& x,
/// Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs.
/// This function also collapses the Jacobian matrices into one like collapsJacs().
inline
AutoDiffBlock<double>
vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
{
typedef AutoDiffBlock<double> ADB;
if (x.empty()) {
return ADB::null();
}
// Count sizes, nonzeros.
const int nx = x.size();
int size = 0;
int nnz = 0;
int elem_with_deriv = -1;
int num_blocks = 0;
for (int elem = 0; elem < nx; ++elem) {
size += x[elem].size();
if (x[elem].derivative().empty()) {
// No nnz contributions from this element.
continue;
} else {
if (elem_with_deriv == -1) {
elem_with_deriv = elem;
num_blocks = x[elem].numBlocks();
}
}
if (x[elem].blockPattern() != x[elem_with_deriv].blockPattern()) {
OPM_THROW(std::runtime_error, "vertcatCollapseJacs(): all arguments must have the same block pattern");
}
for (int block = 0; block < num_blocks; ++block) {
nnz += x[elem].derivative()[block].nonZeros();
}
}
int num_cols = 0;
for (int block = 0; block < num_blocks; ++block) {
num_cols += x[elem_with_deriv].derivative()[block].cols();
}
// Build value for result.
ADB::V val(size);
int pos = 0;
for (int elem = 0; elem < nx; ++elem) {
const int loc_size = x[elem].size();
val.segment(pos, loc_size) = x[elem].value();
pos += loc_size;
}
assert(pos == size);
// Return a constant if we have no derivatives at all.
if (num_blocks == 0) {
return ADB::constant(std::move(val));
}
// Set up for batch insertion of all Jacobian elements.
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> t;
t.reserve(nnz);
int block_row_start = 0;
for (int elem = 0; elem < nx; ++elem) {
int block_col_start = 0;
if (!x[elem].derivative().empty()) {
for (int block = 0; block < num_blocks; ++block) {
const ADB::M& jac = x[elem].derivative()[block];
for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) {
for (ADB::M::InnerIterator i(jac, k); i ; ++i) {
t.push_back(Tri(i.row() + block_row_start,
i.col() + block_col_start,
i.value()));
}
}
block_col_start += jac.cols();
}
}
block_row_start += x[elem].size();
}
// Build final jacobian.
std::vector<ADB::M> jac(1);
jac[0] = Eigen::SparseMatrix<double>(size, num_cols);
jac[0].reserve(nnz);
jac[0].setFromTriplets(t.begin(), t.end());
// Use move semantics to return result efficiently.
return ADB::function(std::move(val), std::move(jac));
}
class Span
{
public:

6
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@@ -463,11 +463,7 @@ namespace Opm
L.setFromTriplets(t.begin(), t.end());
// Combine in single block.
ADB total_residual = std::move(eqs[0]);
for (int phase = 1; phase < num_phases; ++phase) {
total_residual = vertcat(total_residual, eqs[phase]);
}
total_residual = collapseJacs(total_residual);
ADB total_residual = vertcatCollapseJacs(eqs);
// Create output as product of L with equations.
A = L * total_residual.derivative()[0];

146
tests/test_autodiffhelpers.cpp Normal file
View File

@@ -0,0 +1,146 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE AutoDiffHelpersTest
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <boost/test/unit_test.hpp>
using namespace Opm;
namespace {
template <typename Scalar>
bool
operator ==(const Eigen::SparseMatrix<Scalar>& A,
const Eigen::SparseMatrix<Scalar>& B)
{
// Two SparseMatrices are equal if
// 0) They have the same ordering (enforced by equal types)
// 1) They have the same outer and inner dimensions
// 2) They have the same number of non-zero elements
// 3) They have the same sparsity structure
// 4) The non-zero elements are equal
// 1) Outer and inner dimensions
bool eq = (A.outerSize() == B.outerSize());
eq = eq && (A.innerSize() == B.innerSize());
// 2) Equal number of non-zero elements
eq = eq && (A.nonZeros() == B.nonZeros());
for (typename Eigen::SparseMatrix<Scalar>::Index
k0 = 0, kend = A.outerSize(); eq && (k0 < kend); ++k0) {
for (typename Eigen::SparseMatrix<Scalar>::InnerIterator
iA(A, k0), iB(B, k0); eq && (iA && iB); ++iA, ++iB) {
// 3) Sparsity structure
eq = (iA.row() == iB.row()) && (iA.col() == iB.col());
// 4) Equal non-zero elements
eq = eq && (iA.value() == iB.value());
}
}
return eq;
// Note: Investigate implementing this operator as
// return A.cwiseNotEqual(B).count() == 0;
}
}
BOOST_AUTO_TEST_CASE(vertcatCollapseJacsTest)
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
// We will build a system with the block structure
// { 2, 0, 1 } (total of three columns) and { 1, 2, 1 } (row) sizes.
//
// value jacobians
// 10 1 2 | 3
// ----------------------------------
// 11 0 0 | 0 (empty jacobian)
// 12 0 0 | 0
// ----------------------------------
// 13 4 5 | 6
std::vector<ADB> v;
{
// First block.
V val(1);
val << 10;
std::vector<M> jacs(3);
jacs[0] = M(1, 2);
jacs[1] = M(1, 0);
jacs[2] = M(1, 1);
jacs[0].insert(0, 0) = 1.0;
jacs[0].insert(0, 1) = 2.0;
jacs[2].insert(0, 0) = 3.0;
v.push_back(ADB::function(val, jacs));
}
{
// Second block (with empty jacobian).
V val(2);
val << 11, 12;
v.push_back(ADB::constant(val));
}
{
// Third block.
V val(1);
val << 13;
std::vector<M> jacs(3);
jacs[0] = M(1, 2);
jacs[1] = M(1, 0);
jacs[2] = M(1, 1);
jacs[0].insert(0, 0) = 4.0;
jacs[0].insert(0, 1) = 5.0;
jacs[2].insert(0, 0) = 6.0;
v.push_back(ADB::function(val, jacs));
}
std::vector<int> expected_block_pattern{ 2, 0, 1 };
BOOST_CHECK(v[0].blockPattern() == expected_block_pattern);
// Call vertcatCollapseJacs().
const ADB x = vertcatCollapseJacs(v);
// Build expected results.
V expected_val(4);
expected_val << 10, 11, 12, 13;
M expected_jac(4, 3);
expected_jac.insert(0, 0) = 1.0;
expected_jac.insert(0, 1) = 2.0;
expected_jac.insert(0, 2) = 3.0;
expected_jac.insert(3, 0) = 4.0;
expected_jac.insert(3, 1) = 5.0;
expected_jac.insert(3, 2) = 6.0;
// Compare.
BOOST_CHECK((x.value() == expected_val).all());
BOOST_CHECK(x.derivative()[0] == expected_jac);
}