From 35b34f1b3af5e7ed28e5123b9071aaea635e30af Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 8 Mar 2016 10:04:19 +0100 Subject: [PATCH] Add pow() for constant base raised to variable exponent in AutoDiffBlock - associated tests are added - this PR also contains some clean up --- opm/autodiff/AutoDiffBlock.hpp | 144 ++++++++++++++++++--------------- tests/test_block.cpp | 23 ++++-- 2 files changed, 98 insertions(+), 69 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index b0c2a4d8a..71c42b746 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -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 AutoDiffBlock pow(const AutoDiffBlock& base, - const typename AutoDiffBlock::V& exp) + const double exponent) { - const int num_elem = base.value().size(); - typename AutoDiffBlock::V val (num_elem); - typename AutoDiffBlock::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::M derivative_diag(derivative.matrix().asDiagonal()); - - std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); - for (int block = 0; block < base.numBlocks(); block++) { - fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]); - } - - return AutoDiffBlock::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 - AutoDiffBlock pow(const AutoDiffBlock& base, - const double exp) - { - const typename AutoDiffBlock::V val = base.value().pow(exp); - const typename AutoDiffBlock::V derivative = exp * base.value().pow(exp - 1.0); + const typename AutoDiffBlock::V val = base.value().pow(exponent); + const typename AutoDiffBlock::V derivative = exponent * base.value().pow(exponent - 1.0); const typename AutoDiffBlock::M derivative_diag(derivative.matrix().asDiagonal()); std::vector< typename AutoDiffBlock::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 + AutoDiffBlock pow(const AutoDiffBlock& base, + const typename AutoDiffBlock::V& exponent) + { + // Add trivial derivatives and use the AD pow function + return pow(base,AutoDiffBlock::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 + AutoDiffBlock pow(const typename AutoDiffBlock::V& base, + const AutoDiffBlock& exponent) + { + // Add trivial derivatives and use the AD pow function + return pow(AutoDiffBlock::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 + */ + template AutoDiffBlock pow(const AutoDiffBlock& base, - const AutoDiffBlock& exp) - { + const AutoDiffBlock& exponent) + { const int num_elem = base.value().size(); - assert(exp.value().size() == num_elem); + assert(exponent.size() == num_elem); typename AutoDiffBlock::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::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::M > jac1 (base.numBlocks()); - const typename AutoDiffBlock::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::M > jac (num_blocks); + + if ( !exponent.derivative().empty() ) { + typename AutoDiffBlock::V der1 = val; + for (int i = 0; i < num_elem; ++i) { + der1[i] *= std::log(base.value()[i]); + } + std::vector< typename AutoDiffBlock::M > jac1 (exponent.numBlocks()); + const typename AutoDiffBlock::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::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::M > jac2 (base.numBlocks()); - const typename AutoDiffBlock::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::M > jac (base.numBlocks()); - for (int block = 0; block < base.numBlocks(); block++) { - jac[block] = jac1[block] + jac2[block]; + + if ( !base.derivative().empty() ) { + typename AutoDiffBlock::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::M > jac2 (base.numBlocks()); + const typename AutoDiffBlock::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::function(std::move(val), std::move(jac)); diff --git a/tests/test_block.cpp b/tests/test_block.cpp index 48cbb3fce..ad3eda5fc 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -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); + } }