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 pow() for constant base raised to variable exponent in AutoDiffBlock

Browse Source 
- associated tests are added
- this PR also contains some clean up
This commit is contained in:
Tor Harald Sandve 2016-03-08 10:04:19 +01:00
parent b2e02f6d2b
commit 35b34f1b3a
2 changed files with 98 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

144
opm/autodiff/AutoDiffBlock.hpp
View File

@ -624,49 +624,20 @@ namespace Opm
return rhs * lhs; // Commutative operation.
}
/**
* @brief Computes the value of base raised to the power of exp elementwise
* @brief Computes the value of base raised to the power of exponent
*
* @param base The AD forward block
* @param exp array of exponents
* @return The value of base raised to the power of exp elementwise
* @param exponent double
* @return The value of base raised to the power of exponent
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const typename AutoDiffBlock<Scalar>::V& exp)
const double exponent)
{
const int num_elem = base.value().size();
typename AutoDiffBlock<Scalar>::V val (num_elem);
typename AutoDiffBlock<Scalar>::V derivative = exp;
assert(exp.size() == num_elem);
for (int i = 0; i < num_elem; ++i) {
val[i] = std::pow(base.value()[i], exp[i]);
derivative[i] *= std::pow(base.value()[i], exp[i] - 1.0);
}
const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
for (int block = 0; block < base.numBlocks(); block++) {
fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]);
}
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 AD forward block
* @param exp exponent
* @return The value of base raised to the power of exp
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const double exp)
{
const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exp);
const typename AutoDiffBlock<Scalar>::V derivative = exp * base.value().pow(exp - 1.0);
const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exponent);
const typename AutoDiffBlock<Scalar>::V derivative = exponent * base.value().pow(exponent - 1.0);
const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
@ -678,44 +649,91 @@ namespace Opm
}
/**
* @brief Computes the value of base raised to the power of exp
* @brief Computes the value of base raised to the power of exponent
*
* @param base The AD forward block
* @param exponent Array of exponents
* @return The value of base raised to the power of exponent elementwise
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const typename AutoDiffBlock<Scalar>::V& exponent)
{
// Add trivial derivatives and use the AD pow function
return pow(base,AutoDiffBlock<Scalar>::constant(exponent));
}
/**
* @brief Computes the value of base raised to the power of exponent
*
* @param base Array of base values
* @param exponent The AD forward block
* @return The value of base raised to the power of exponent elementwise
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const typename AutoDiffBlock<Scalar>::V& base,
const AutoDiffBlock<Scalar>& exponent)
{
// Add trivial derivatives and use the AD pow function
return pow(AutoDiffBlock<Scalar>::constant(base),exponent);
}
/**
* @brief Computes the value of base raised to the power of exponent
*
* @param base The base AD forward block
* @param exp The exponent AD forward block
* @param exponent The exponent AD forward block
* @return The value of base raised to the power of exp
*/ template <typename Scalar>
*/
template <typename Scalar>
AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
const AutoDiffBlock<Scalar>& exp)
{
const AutoDiffBlock<Scalar>& exponent)
{
const int num_elem = base.value().size();
assert(exp.value().size() == num_elem);
assert(exponent.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]);
val[i] = std::pow(base.value()[i], exponent.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]);
// (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' = der1 + der2
// if f' is empty only der1 is calculated
// if g' is empty only der2 is calculated
// if f' and g' are non empty they should have the same size
int num_blocks = std::max (base.numBlocks(), exponent.numBlocks());
if (!base.derivative().empty() && !exponent.derivative().empty()) {
assert(exponent.numBlocks() == base.numBlocks());
}
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]);
std::vector< typename AutoDiffBlock<Scalar>::M > jac (num_blocks);
if ( !exponent.derivative().empty() ) {
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 (exponent.numBlocks());
const typename AutoDiffBlock<Scalar>::M der1_diag(der1.matrix().asDiagonal());
for (int block = 0; block < exponent.numBlocks(); block++) {
fastSparseProduct(der1_diag, exponent.derivative()[block], jac1[block]);
jac[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];
if ( !base.derivative().empty() ) {
typename AutoDiffBlock<Scalar>::V der2 = exponent.value();
for (int i = 0; i < num_elem; ++i) {
der2[i] *= std::pow(base.value()[i], exponent.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]);
if (!exponent.derivative().empty()) {
jac[block] += jac2[block];
} else {
jac[block] = jac2[block];
}
}
}
return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));

23
tests/test_block.cpp
View File

@ -340,20 +340,31 @@ BOOST_AUTO_TEST_CASE(Pow)
ADB compare = pick1 * xx + pick2 * xxx + pick3 * xpowhalf;
checkClose(xpowyval, compare, tolerance);
// test exp = ADB
ADB xpowy = Opm::pow(x,y);
// test exponent = ADB::V and base = ADB
ADB xvalpowy = Opm::pow(x.value(),y);
// the value and the first jacobian should be equal to the xpowyval
// the value should be equal to xpowyval
// the first jacobian should be trivial
// 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);
BOOST_CHECK_CLOSE(xvalpowy.value()[i], xpowyval.value()[i], tolerance);
BOOST_CHECK_CLOSE(xvalpowy.derivative()[0].coeff(i,i), 0.0, tolerance);
BOOST_CHECK_CLOSE(xvalpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4);
}
// test exp = ADB
ADB xpowy = Opm::pow(x,y);
// the first jacobian should be equal to the xpowyval
// the second jacobian should be equal to the xvalpowy
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), xvalpowy.derivative()[1].coeff(i,i), tolerance);
}
}