From b2e02f6d2b465e734bec7300e77277cbc4b1bede Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Mon, 7 Mar 2016 14:23:45 +0100
Subject: [PATCH] Add power method for general f^g in the AutoDiffBlock

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
---
 opm/autodiff/AutoDiffBlock.hpp | 45 ++++++++++++++++++++++++++++++++++
 tests/test_block.cpp           | 35 +++++++++++++++++++++++++-
 2 files changed, 79 insertions(+), 1 deletion(-)

diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp
index 66901ea4e..b0c2a4d8a 100644
--- a/opm/autodiff/AutoDiffBlock.hpp
+++ b/opm/autodiff/AutoDiffBlock.hpp
@@ -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
diff --git a/tests/test_block.cpp b/tests/test_block.cpp
index 30928d400..48cbb3fce 100644
--- a/tests/test_block.cpp
+++ b/tests/test_block.cpp
@@ -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);
+    }
+
+
 }