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

Add power method for general f^g in the AutoDiffBlock

Browse Source 
A power method where both f and g are ADB variables is added
using the general derivative rule
(f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f'

Tests are added to test_block.cpp
This commit is contained in:
Tor Harald Sandve 2016-03-07 14:23:45 +01:00
parent 96938548c1
commit b2e02f6d2b
2 changed files with 79 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
opm/autodiff/AutoDiffBlock.hpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2016 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -676,6 +677,50 @@ namespace Opm
return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
}
/**
* @brief Computes the value of base raised to the power of exp
*
* @param base The base AD forward block
* @param exp The exponent AD forward block
* @return The value of base raised to the power of exp
*/ template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const AutoDiffBlock<Scalar>& exp)
{
const int num_elem = base.value().size();
assert(exp.value().size() == num_elem);
typename AutoDiffBlock<Scalar>::V val (num_elem);
for (int i = 0; i < num_elem; ++i) {
val[i] = std::pow(base.value()[i], exp.value()[i]);
}
// (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f'
typename AutoDiffBlock<Scalar>::V der1 = val;
for (int i = 0; i < num_elem; ++i) {
der1[i] *= std::log(base.value()[i]);
}
std::vector< typename AutoDiffBlock<Scalar>::M > jac1 (base.numBlocks());
const typename AutoDiffBlock<Scalar>::M der1_diag(der1.matrix().asDiagonal());
for (int block = 0; block < base.numBlocks(); block++) {
fastSparseProduct(der1_diag, exp.derivative()[block], jac1[block]);
}
typename AutoDiffBlock<Scalar>::V der2 = exp.value();
for (int i = 0; i < num_elem; ++i) {
der2[i] *= std::pow(base.value()[i], exp.value()[i] - 1.0);
}
std::vector< typename AutoDiffBlock<Scalar>::M > jac2 (base.numBlocks());
const typename AutoDiffBlock<Scalar>::M der2_diag(der2.matrix().asDiagonal());
for (int block = 0; block < base.numBlocks(); block++) {
fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]);
}
std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
for (int block = 0; block < base.numBlocks(); block++) {
jac[block] = jac1[block] + jac2[block];
}
return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
}
} // namespace Opm

35
tests/test_block.cpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2016 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -293,7 +294,7 @@ BOOST_AUTO_TEST_CASE(Pow)
vx << 0.2, 1.2, 13.4;
ADB::V vy(3);
vy << 1.0, 2.2, 3.4;
vy << 2.0, 3.0, 0.5;
std::vector<ADB::V> vals{ vx, vy };
std::vector<ADB> vars = ADB::variables(vals);
@ -303,6 +304,7 @@ BOOST_AUTO_TEST_CASE(Pow)
const double tolerance = 1e-14;
// test exp = double
const ADB xx = x * x;
ADB xxpow2 = Opm::pow(x,2.0);
checkClose(xxpow2, xx, tolerance);
@ -322,6 +324,37 @@ BOOST_AUTO_TEST_CASE(Pow)
for (int i = 0 ; i < 3; ++i){
BOOST_CHECK_CLOSE(xpowhalf.value()[i], x_sqrt[i], 1e-4);
}
// test exp = ADB::V
ADB xpowyval = Opm::pow(x,y.value());
// each of the component of y is tested in the test above
// we compare with the results from the above tests.
ADB::V pick1(3);
pick1 << 1,0,0;
ADB::V pick2(3);
pick2 << 0,1,0;
ADB::V pick3(3);
pick3 << 0,0,1;
ADB compare = pick1 * xx + pick2 * xxx + pick3 * xpowhalf;
checkClose(xpowyval, compare, tolerance);
// test exp = ADB
ADB xpowy = Opm::pow(x,y);
// the value and the first jacobian should be equal to the xpowyval
// the second jacobian is hand calculated
// log(0.2)*0.2^2.0, log(1.2) * 1.2^3.0, log(13.4) * 13.4^0.5
ADB::V jac2(3);
jac2 << -0.0643775165 , 0.315051650 , 9.50019208855;
for (int i = 0 ; i < 3; ++i){
BOOST_CHECK_CLOSE(xpowy.value()[i], xpowyval.value()[i], tolerance);
BOOST_CHECK_CLOSE(xpowy.derivative()[0].coeff(i,i), xpowyval.derivative()[0].coeff(i,i), tolerance);
BOOST_CHECK_CLOSE(xpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4);
}
}