Merge remote-tracking branch 'atgeirr/master'

Conflicts:
	sim_simple.cpp

AutoDiffBlock.hpp

@@ -156,6 +156,17 @@ namespace AutoDiff
             return jac_.size();
         }
 
+        /// Sizes (number of columns) of Jacobian blocks.
+        std::vector<int> blockPattern() const
+        {
+            const int nb = numBlocks();
+            std::vector<int> bp(nb);
+            for (int block = 0; block < nb; ++block) {
+                bp[block] = jac_[block].cols();
+            }
+            return bp;
+        }
+
         /// Function value
         const V& value() const
         {
@@ -194,6 +205,7 @@ namespace AutoDiff
         return fw.print(os);
     }
 
+
     /// Multiply with sparse matrix from the left.
     template <typename Scalar>
     ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::M& lhs,
@@ -210,6 +222,70 @@ namespace AutoDiff
     }
 
 
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::V& lhs,
+                                   const ForwardBlock<Scalar>& rhs)
+    {
+        return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator*(const ForwardBlock<Scalar>& lhs,
+                                   const typename ForwardBlock<Scalar>::V& rhs)
+    {
+        return rhs * lhs; // Commutative operation.
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator+(const typename ForwardBlock<Scalar>::V& lhs,
+                                   const ForwardBlock<Scalar>& rhs)
+    {
+        return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator+(const ForwardBlock<Scalar>& lhs,
+                                   const typename ForwardBlock<Scalar>::V& rhs)
+    {
+        return rhs + lhs; // Commutative operation.
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator-(const typename ForwardBlock<Scalar>::V& lhs,
+                                   const ForwardBlock<Scalar>& rhs)
+    {
+        return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator-(const ForwardBlock<Scalar>& lhs,
+                                   const typename ForwardBlock<Scalar>::V& rhs)
+    {
+        return lhs - ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator/(const typename ForwardBlock<Scalar>::V& lhs,
+                                   const ForwardBlock<Scalar>& rhs)
+    {
+        return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
+    }
+
+
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator/(const ForwardBlock<Scalar>& lhs,
+                                   const typename ForwardBlock<Scalar>::V& rhs)
+    {
+        return lhs / ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
+    }
+
+
 } // namespace Autodiff
 
 

AutoDiffBlockArma.hpp

@@ -0,0 +1,216 @@
+/*
+  Copyright 2013 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/>.
+*/
+
+#ifndef OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED
+#define OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED
+
+
+// #include "AutoDiff.hpp"
+// #include <Eigen/Eigen>
+// #include <Eigen/Sparse>
+#include <armadillo>
+#include <vector>
+#include <cassert>
+
+namespace AutoDiff
+{
+
+    template <typename Scalar>
+    class ForwardBlock
+    {
+    public:
+        /// Underlying types for scalar vectors and jacobians.
+        typedef arma::Col<Scalar> V;
+        typedef arma::SpMat<Scalar> M;
+
+        /// Named constructor pattern used here.
+        static ForwardBlock constant(const V& val, const std::vector<int>& blocksizes)
+        {
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // For constants, all jacobian blocks are zero.
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            return ForwardBlock(val, jac);
+        }
+
+        static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
+        {
+            std::vector<M> jac;
+            const int num_elem = val.size();
+            const int num_blocks = blocksizes.size();
+            // First, set all jacobian blocks to zero...
+            jac.resize(num_blocks);
+            for (int i = 0; i < num_blocks; ++i) {
+                jac[i] = M(num_elem, blocksizes[i]);
+            }
+            // ... then set the one corrresponding to this variable to identity.
+            assert(blocksizes[index] == num_elem);
+            jac[index].eye(num_elem, num_elem);
+            return ForwardBlock(val, jac);
+        }
+
+        static ForwardBlock function(const V& val, const std::vector<M>& jac)
+        {
+            return ForwardBlock(val, jac);
+        }
+
+        /// Operator +
+        ForwardBlock operator+(const ForwardBlock& rhs)
+        {
+            std::vector<M> jac = jac_;
+            assert(numBlocks() == rhs.numBlocks());
+            int num_blocks = numBlocks();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac[block].n_rows == rhs.jac_[block].n_rows);
+                assert(jac[block].n_cols == rhs.jac_[block].n_cols);
+                jac[block] += rhs.jac_[block];
+            }
+            return function(val_ + rhs.val_, jac);
+        }
+
+        /// Operator -
+        ForwardBlock operator-(const ForwardBlock& rhs)
+        {
+            std::vector<M> jac = jac_;
+            assert(numBlocks() == rhs.numBlocks());
+            int num_blocks = numBlocks();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac[block].n_rows == rhs.jac_[block].n_rows);
+                assert(jac[block].n_cols == rhs.jac_[block].n_cols);
+                jac[block] -= rhs.jac_[block];
+            }
+            return function(val_ - rhs.val_, jac);
+        }
+
+        /// Operator *
+        ForwardBlock operator*(const ForwardBlock& rhs)
+        {
+            int num_blocks = numBlocks();
+            std::vector<M> jac(num_blocks);
+            assert(numBlocks() == rhs.numBlocks());
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac_[block].n_rows == rhs.jac_[block].n_rows);
+                assert(jac_[block].n_cols == rhs.jac_[block].n_cols);
+                jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
+            }
+            return function(val_ * rhs.val_, jac);
+        }
+
+        /// Operator /
+        ForwardBlock operator/(const ForwardBlock& rhs)
+        {
+            int num_blocks = numBlocks();
+            std::vector<M> jac(num_blocks);
+            assert(numBlocks() == rhs.numBlocks());
+            typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
+            D D1 = val_.matrix().asDiagonal();
+            D D2 = rhs.val_.matrix().asDiagonal();
+            D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(jac_[block].n_rows == rhs.jac_[block].n_rows);
+                assert(jac_[block].n_cols == rhs.jac_[block].n_cols);
+                jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
+            }
+            return function(val_ / rhs.val_, jac);
+        }
+        /// I/O.
+        template <class Ostream>
+        Ostream&
+        print(Ostream& os) const
+        {
+            int num_blocks = jac_.size();
+            os << "Value =\n" << val_ << "\n\nJacobian =\n";
+            for (int i = 0; i < num_blocks; ++i) {
+                os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
+            }
+            return os;
+        }
+
+        /// Number of variables or Jacobian blocks.
+        int numBlocks() const
+        {
+            return jac_.size();
+        }
+
+        /// Function value
+        const V& value() const
+        {
+            return val_;
+        }
+
+        /// Function derivatives
+        const std::vector<M>& derivative() const
+        {
+            return jac_;
+        }
+
+    private:
+        ForwardBlock(const V& val,
+                     const std::vector<M>& jac)
+            : val_(val), jac_(jac)
+        {
+#ifndef NDEBUG
+            const int num_elem = val_.size();
+            const int num_blocks = jac_.size();
+            for (int block = 0; block < num_blocks; ++block) {
+                assert(num_elem == jac_[block].n_rows);
+            }
+#endif
+        }
+
+        V val_;
+        std::vector<M> jac_;
+    };
+
+
+    template <class Ostream, typename Scalar>
+    Ostream&
+    operator<<(Ostream& os, const ForwardBlock<Scalar>& fw)
+    {
+        return fw.print(os);
+    }
+
+    /// Multiply with sparse matrix from the left.
+    template <typename Scalar>
+    ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::M& lhs,
+                                   const ForwardBlock<Scalar>& rhs)
+    {
+        int num_blocks = rhs.numBlocks();
+        std::vector<typename ForwardBlock<Scalar>::M> jac(num_blocks);
+        assert(lhs.n_cols == rhs.value().n_rows);
+        for (int block = 0; block < num_blocks; ++block) {
+            jac[block] = lhs*rhs.derivative()[block];
+        }
+        typename ForwardBlock<Scalar>::V val = lhs*rhs.value().matrix();
+        return ForwardBlock<Scalar>::function(val, jac);
+    }
+
+
+} // namespace Autodiff
+
+
+
+#endif // OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED

sim_simple.cpp

@@ -55,8 +55,8 @@ struct HelperOps
     typedef AutoDiff::ForwardBlock<double>::V V;
 
     /// A list of internal faces.
-    typedef Eigen::Array<int, Eigen::Dynamic, 1> iFaces;
-    iFaces internal_faces;
+    typedef Eigen::Array<int, Eigen::Dynamic, 1> IFaces;
+    IFaces internal_faces;
 
     /// Extract for each face the difference of its adjacent cells'values.
     M ngrad;
@@ -172,13 +172,13 @@ namespace {
         select(const typename ADB::V&  press,
                const std::vector<ADB>& xc   ) const
         {
-            typedef HelperOps::iFaces::Index ifIndex;
-            const ifIndex nif = h_.internal_faces.size();
+            typedef HelperOps::IFaces::Index IFIndex;
+            const IFIndex nif = h_.internal_faces.size();
 
             // Define selector structure.
             typedef typename Eigen::Triplet<Scalar> Triplet;
             std::vector<Triplet> s;  s.reserve(nif);
-            for (ifIndex i = 0; i < nif; ++i) {
+            for (IFIndex i = 0; i < nif; ++i) {
                 const int f  = h_.internal_faces[i];
                 const int c1 = g_.face_cells[2*f + 0];
                 const int c2 = g_.face_cells[2*f + 1];
@@ -276,27 +276,27 @@ int main()
 
     Opm::time::StopWatch clock;
     clock.start();
-    Opm::GridManager gm(3,3);//(50, 50, 10);
+    const Opm::GridManager gm(3,3);//(50, 50, 10);
     const UnstructuredGrid& grid = *gm.c_grid();
     using namespace Opm::unit;
     using namespace Opm::prefix;
-    // Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
-    //                                  { 1000.0, 800.0 },
-    //                                  { 1.0*centi*Poise, 5.0*centi*Poise },
-    //                                  0.2, 100*milli*darcy,
-    //                                  grid.dimensions, grid.number_of_cells);
-    // Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
-    //                                  { 1000.0, 1000.0 },
-    //                                  { 1.0, 1.0 },
-    //                                  1.0, 1.0,
-    //                                  grid.dimensions, grid.number_of_cells);
-    Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
-                                     { 1000.0, 1000.0 },
-                                     { 1.0, 30.0 },
-                                     1.0, 1.0,
-                                     grid.dimensions, grid.number_of_cells);
+    // const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
+    //                                        { 1000.0, 800.0 },
+    //                                        { 1.0*centi*Poise, 5.0*centi*Poise },
+    //                                        0.2, 100*milli*darcy,
+    //                                        grid.dimensions, grid.number_of_cells);
+    // const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
+    //                                        { 1000.0, 1000.0 },
+    //                                        { 1.0, 1.0 },
+    //                                        1.0, 1.0,
+    //                                        grid.dimensions, grid.number_of_cells);
+    const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear,
+                                           { 1000.0, 1000.0 },
+                                           { 1.0, 30.0 },
+                                           1.0, 1.0,
+                                           grid.dimensions, grid.number_of_cells);
     std::vector<double> htrans(grid.cell_facepos[grid.number_of_cells]);
-    tpfa_htrans_compute((UnstructuredGrid*)&grid, props.permeability(), htrans.data());
+    tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
     // std::vector<double> trans(grid.number_of_faces);
     V trans_all(grid.number_of_faces);
     tpfa_trans_compute((UnstructuredGrid*)&grid, htrans.data(), trans_all.data());
@@ -308,7 +308,7 @@ int main()
     std::cerr << "Opm core " << clock.secsSinceLast() << std::endl;
 
     // Define neighbourhood-derived operator matrices.
-    HelperOps ops(grid);
+    const HelperOps ops(grid);
     const int num_internal = ops.internal_faces.size();
     V transi(num_internal);
     for (int fi = 0; fi < num_internal; ++fi) {
@@ -334,15 +334,15 @@ int main()
     // totmob - explicit as well
     TwoCol kr(nc, 2);
     props.relperm(nc, s0.data(), allcells.data(), kr.data(), 0);
-    V krw = kr.leftCols<1>();
-    V kro = kr.rightCols<1>();
+    const V krw = kr.leftCols<1>();
+    const V kro = kr.rightCols<1>();
     const double* mu = props.viscosity();
-    V totmob = krw/mu[0] + kro/mu[1];
-    V totmobf = (ops.caver*totmob.matrix()).array();
+    const V totmob = krw/mu[0] + kro/mu[1];
+    const V totmobf = (ops.caver*totmob.matrix()).array();
 
     // Mobility-weighted transmissibilities per internal face.
     // Still explicit, and no upwinding!
-    V mobtransf = totmobf*transi;
+    const V mobtransf = totmobf*transi;
 
     std::cerr << "Property arrays " << clock.secsSinceLast() << std::endl;
 
@@ -351,18 +351,14 @@ int main()
     p0.fill(200*Opm::unit::barsa);
 
     // First actual AD usage: defining pressure variable.
-    std::vector<int> block_pattern = { nc };
-    // Could actually write { nc } instead of block_pattern below,
+    const std::vector<int> bpat = { nc };
+    // Could actually write { nc } instead of bpat below,
     // but we prefer a named variable since we will repeat it.
-    ADB p = ADB::variable(0, p0, block_pattern);
-    ADB ngradp = ops.ngrad*p;
+    const ADB p = ADB::variable(0, p0, bpat);
+    const ADB ngradp = ops.ngrad*p;
     // We want flux = totmob*trans*(p_i - p_j) for the ij-face.
-    // We only need to multiply mobtransf and pdiff_face,
-    // but currently multiplication with constants is not in,
-    // so we define an AD constant to multiply with.
-    ADB mobtransf_ad = ADB::constant(mobtransf, block_pattern);
-    ADB flux = mobtransf_ad*ngradp;
-    ADB residual = ops.div*flux - ADB::constant(q, block_pattern);
+    const ADB flux = mobtransf*ngradp;
+    const ADB residual = ops.div*flux - q;
     std::cerr << "Construct AD residual " << clock.secsSinceLast() << std::endl;
 
     // It's the residual we want to be zero. We know it's linear in p,
@@ -374,20 +370,21 @@ int main()
     // residual.derived()[0].
 
     Eigen::UmfPackLU<M> solver;
-    M matr = residual.derivative()[0];
-    matr.coeffRef(0,0) *= 2.0;
-    matr.makeCompressed();
-    solver.compute(matr);
+    M pmatr = residual.derivative()[0];
+    pmatr.coeffRef(0,0) *= 2.0;
+    pmatr.makeCompressed();
+    solver.compute(pmatr);
     if (solver.info() != Eigen::Success) {
         std::cerr << "Pressure/flow Jacobian decomposition error\n";
         return EXIT_FAILURE;
     }
-    Eigen::VectorXd x = solver.solve(residual.value().matrix());
+    // const Eigen::VectorXd dp = solver.solve(residual.value().matrix());
+    const V dp = solver.solve(residual.value().matrix()).array();
     if (solver.info() != Eigen::Success) {
         std::cerr << "Pressure/flow solve failure\n";
         return EXIT_FAILURE;
     }
-    V p1 = p0 - x.array();
+    const V p1 = p0 - dp;
     std::cerr << "Solve " << clock.secsSinceLast() << std::endl;
     // std::cout << p1 << std::endl;
 
@@ -401,69 +398,52 @@ int main()
     // and f_w is (for now) based on averaged mobilities, not upwind.
 
     double res_norm = 1e100;
-    V s1 =  /*s0.leftCols<1>()*/0.5*V::Ones(nc,1); // Initial guess.
-    UpwindSelector<double> upws(grid, ops);
-    const ADB nkdp = (ADB::constant(transi                 , block_pattern) *
-                      ADB::constant(ops.ngrad * p1.matrix(), block_pattern));
-    const ADB s00 = ADB::constant(s0.leftCols<1>(), block_pattern);
-    const std::vector<ADB> pmobc0 = phaseMobility<ADB>(props, allcells, s00.value());
-    const std::vector<ADB> pmobf0 = upws.select(p1, pmobc0);
-    const std::vector<ADB::M> null = { ADB::M(transi.size(), nc) };
-    const ADB dflux = (ADB::function((pmobf0[0] + pmobf0[1]).value(), null) *
-                       ADB::function(nkdp.value()                   , null));
+    const V sw0 = s0.leftCols<1>();
+    // V sw1 = sw0;
+    V sw1 = 0.5*V::Ones(nc,1);
+    const UpwindSelector<double> upws(grid, ops);
+    const V nkdp = transi * (ops.ngrad * p1.matrix()).array();
+    const V dflux = totmobf * nkdp;
+    const V pv = Eigen::Map<const V>(props.porosity(), nc, 1)
+        * Eigen::Map<const V>(grid.cell_volumes, nc, 1);
+    const double dt = 0.0005;
+    const V dtpv = dt/pv;
+    const V qneg = q.min(V::Zero(nc,1));
+    const V qpos = q.max(V::Zero(nc,1));
 
     std::cout.setf(std::ios::scientific);
     std::cout.precision(16);
 
     int it = 0;
     do {
-        const std::vector<int>& bp = block_pattern;
-        ADB s = ADB::variable(0, s1, bp);
-        const double dt = 0.0005;
-        V pv = Eigen::Map<const V>(props.porosity(), nc, 1)
-            * Eigen::Map<const V>(grid.cell_volumes, nc, 1);
-        V dtpv = dt/pv;
-        // std::cout << dtpv;
-        std::vector<ADB> pmobc = phaseMobility<ADB>(props, allcells, s.value());
-        std::vector<ADB> pmobf = upws.select(p1, pmobc);
-
-        ADB fw_cell = fluxFunc(pmobc);
+        const ADB sw = ADB::variable(0, sw1, bpat);
+        const std::vector<ADB> pmobc = phaseMobility<ADB>(props, allcells, sw.value());
+        const std::vector<ADB> pmobf = upws.select(p1, pmobc);
+        const ADB fw_cell = fluxFunc(pmobc);
         const ADB fw_face = fluxFunc(pmobf);
-        ADB flux1 = fw_face * dflux;
-        // std::cout << flux1;
-        V qneg = dtpv*q;
-        V qpos = dtpv*q;
-        // Cheating a bit...
-        qneg[0] = 0.0;
-        qpos[nc-1] = 0.0;
-        ADB qtr_ad = ADB::constant(qpos, bp) + fw_cell*ADB::constant(qneg, bp);
-        ADB transport_residual = s - ADB::constant(s0.leftCols<1>(), bp)
-            + ADB::constant(dtpv, bp)*(ops.div*flux1)
-            - qtr_ad;
+        const ADB flux1 = fw_face * dflux;
+        const ADB qtr_ad = qpos + fw_cell*qneg;
+        const ADB transport_residual = sw - sw0 + dtpv*(ops.div*flux1 - qtr_ad);
         res_norm = transport_residual.value().matrix().norm();
         std::cout << "res_norm[" << it << "] = "
                   << res_norm << std::endl;
 
-        matr = transport_residual.derivative()[0];
-        matr.makeCompressed();
-        // std::cout << transport_residual;
-        solver.compute(matr);
+        M smatr = transport_residual.derivative()[0];
+        smatr.makeCompressed();
+        solver.compute(smatr);
         if (solver.info() != Eigen::Success) {
             std::cerr << "Transport Jacobian decomposition error\n";
             return EXIT_FAILURE;
         }
-        x = solver.solve(transport_residual.value().matrix());
+        const V ds = solver.solve(transport_residual.value().matrix()).array();
         if (solver.info() != Eigen::Success) {
             std::cerr << "Transport solve failure\n";
             return EXIT_FAILURE;
         }
-        // std::cout << x << std::endl;
-        s1 = s.value() - x.array();
+        sw1 = sw.value() - ds;
         std::cerr << "Solve for s[" << it << "]: "
                   << clock.secsSinceLast() << '\n';
-        for (int c = 0; c < nc; ++c) {
-            s1[c] = std::min(1.0, std::max(0.0, s1[c]));
-        }
+        sw1 = sw1.min(V::Ones(nc,1)).max(V::Zero(nc,1));
 
         it += 1;
     } while (res_norm > 1e-7);