Merge remote-tracking branch 'bska/master'

AutoDiffBlock.hpp

@@ -76,7 +76,7 @@ namespace AutoDiff
         }
 
         /// Operator +
-        ForwardBlock operator+(const ForwardBlock& rhs)
+        ForwardBlock operator+(const ForwardBlock& rhs) const
         {
             std::vector<M> jac = jac_;
             assert(numBlocks() == rhs.numBlocks());
@@ -90,7 +90,7 @@ namespace AutoDiff
         }
 
         /// Operator -
-        ForwardBlock operator-(const ForwardBlock& rhs)
+        ForwardBlock operator-(const ForwardBlock& rhs) const
         {
             std::vector<M> jac = jac_;
             assert(numBlocks() == rhs.numBlocks());
@@ -104,7 +104,7 @@ namespace AutoDiff
         }
 
         /// Operator *
-        ForwardBlock operator*(const ForwardBlock& rhs)
+        ForwardBlock operator*(const ForwardBlock& rhs) const
         {
             int num_blocks = numBlocks();
             std::vector<M> jac(num_blocks);
@@ -121,7 +121,7 @@ namespace AutoDiff
         }
 
         /// Operator /
-        ForwardBlock operator/(const ForwardBlock& rhs)
+        ForwardBlock operator/(const ForwardBlock& rhs) const
         {
             int num_blocks = numBlocks();
             std::vector<M> jac(num_blocks);

mrst/sim_simple.m

@@ -0,0 +1,27 @@
+function x = sim_simple(cartDims, physDims, tf, verb)
+   g = computeGeometry(cartGrid(cartDims, physDims));
+
+   rock = struct('perm', repmat(1, [g.cells.num, 1]), ...
+                 'poro', repmat(1, [g.cells.num, 1]));
+
+   T = computeTrans(g, rock);
+
+   fluid = initSimpleFluid('n'  , [   1,   1], ...
+                           'mu' , [   1,  30], ...
+                           'rho', [1000, 800]);
+
+   gravity reset off
+
+   src = addSource([], [1, g.cells.num], [1, -1], 'sat', [ 1, 0 ; 0, 1]);
+
+   s0    = [0, 1];
+   state = incompTPFA(initState(g, [], 0, s0), ...
+                      g, T, fluid, 'src', src);
+
+   if nargin < 4, verb = false; end
+   state = implicitTransport(state, g, tf, rock, fluid, ...
+                             'src', src, 'verbose', verb);
+
+   x = struct('g', g, 'rock', rock, 'T', T, 'fluid', fluid, ...
+              'src', src, 'state', state);
+end

sim_simple.cpp

@@ -55,7 +55,8 @@ struct HelperOps
     typedef AutoDiff::ForwardBlock<double>::V V;
 
     /// A list of internal faces.
-    Eigen::Array<int, Eigen::Dynamic, 1> 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;
@@ -111,11 +112,115 @@ struct HelperOps
     }
 };
 
-/// Returns fw(sw).
+#if !defined(NDEBUG)
+#include <cstdio>
+#endif  // !defined(NDEBUG)
+
+namespace {
+#if !defined(NDEBUG)
+    void
+    printSparseMatrix(const Eigen::SparseMatrix<double>& A,
+                      std::FILE*                         fp)
+    {
+        typedef Eigen::SparseMatrix<double>::Index Index;
+
+        const Index*  const p = A.outerIndexPtr();
+        const Index*  const i = A.innerIndexPtr();
+        const double* const x = A.valuePtr();
+
+        const Index cols = A.outerSize();
+        assert (A.innerSize() == cols);
+
+        for (Index j = 0; j < cols; j++) {
+            for (Index k = p[j]; k < p[j + 1]; k++) {
+                std::fprintf(fp, "%lu %lu %26.18e\n",
+                             static_cast<unsigned long>(i[k] + 1),
+                             static_cast<unsigned long>(j    + 1), x[k]);
+            }
+        }
+    }
+
+    void
+    printSparseMatrix(const Eigen::SparseMatrix<double>& A ,
+                      const char* const                  fn)
+    {
+        std::FILE* fp;
+
+        fp = std::fopen(fn, "w");
+        if (fp != 0) {
+            printSparseMatrix(A, fp);
+        }
+
+        std::fclose(fp);
+    }
+#endif  // !defined(NDEBUG)
+
+    template <typename Scalar>
+    class UpwindSelector {
+    public:
+        typedef AutoDiff::ForwardBlock<Scalar> ADB;
+
+        UpwindSelector(const UnstructuredGrid& g,
+                       const HelperOps&        h)
+            : g_(g), h_(h)
+        {
+        }
+
+        // Upwind selection in absence of counter-current flow (i.e.,
+        // without effects of gravity and/or capillary pressure).
+        std::vector<ADB>
+        select(const typename ADB::V&  press,
+               const std::vector<ADB>& xc   ) const
+        {
+            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) {
+                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];
+
+                assert ((c1 >= 0) && (c2 >= 0));
+
+                const Scalar dp = press[c1] - press[c2];
+                const int c = (! (dp < Scalar(0))) ? c1 : c2;
+
+                s.push_back(Triplet(i, c, Scalar(1)));
+            }
+
+            // Assemble explicit selector operator.
+            typename ADB::M S(nif, g_.number_of_cells);
+            S.setFromTriplets(s.begin(), s.end());
+
+            // Apply selector.
+            //
+            // Absence of counter-current flow means that the same
+            // selector applies to all quantities, 'x', defined per
+            // cell.
+            std::vector<ADB> xf;  xf.reserve(xc.size());
+            for (typename std::vector<ADB>::const_iterator
+                     b = xc.begin(), e = xc.end(); b != e; ++b)
+            {
+                xf.push_back(S * (*b));
+            }
+
+            return xf;
+        }
+
+    private:
+        const UnstructuredGrid& g_;
+        const HelperOps&        h_;
+    };
+}
+
 template <class ADB>
-ADB fluxFunc(const Opm::IncompPropertiesInterface& props,
+std::vector<ADB>
-             const std::vector<int>& cells,
+phaseMobility(const Opm::IncompPropertiesInterface& props,
-             const typename ADB::V& sw)
+              const std::vector<int>& cells,
+              const typename ADB::V& sw)
 {
     typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
     typedef Eigen::Array<double, Eigen::Dynamic, 4, Eigen::RowMajor> FourCol;
@@ -144,11 +249,22 @@ ADB fluxFunc(const Opm::IncompPropertiesInterface& props,
     const double* mu = props.viscosity();
     std::vector<M> dmw = { krwjac/mu[0] };
     std::vector<M> dmo = { krojac/mu[1] };
-    ADB mw_ad = ADB::function(krw/mu[0], dmw);
+
-    ADB mo_ad = ADB::function(kro/mu[1], dmo);
+    std::vector<ADB> pmobc = { ADB::function(krw / mu[0], dmw) ,
-    ADB fw = mw_ad / (mw_ad + mo_ad);
+                               ADB::function(kro / mu[1], dmo) };
-    // std::cout << mw_ad << mo_ad << (mw_ad + mo_ad) << fw;
+    return pmobc;
-    return fw;
+}
+
+/// Returns fw(sw).
+template <class ADB>
+ADB
+fluxFunc(const std::vector<ADB>& m)
+{
+    assert (m.size() == 2);
+
+    ADB f = m[0] / (m[0] + m[1]);
+
+    return f;
 }
 
 
@@ -286,6 +402,20 @@ int main()
 
     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));
+
+    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);
@@ -294,13 +424,12 @@ int main()
             * Eigen::Map<const V>(grid.cell_volumes, nc, 1);
         V dtpv = dt/pv;
         // std::cout << dtpv;
-        V ngradp1 = ops.ngrad*p1.matrix();
+        std::vector<ADB> pmobc = phaseMobility<ADB>(props, allcells, s.value());
-        // std::cout << ngradp1 << std::endl;
+        std::vector<ADB> pmobf = upws.select(p1, pmobc);
-        ADB fw_cell = fluxFunc<ADB>(props, allcells, s.value());
+
-        // std::cout << fw_cell;
+        ADB fw_cell = fluxFunc(pmobc);
-        ADB fw_face = ops.caver*fw_cell; 
+        const ADB fw_face = fluxFunc(pmobf);
-        // std::cout << fw_face;
+        ADB flux1 = fw_face * dflux;
-        ADB flux1 = fw_face*ADB::constant(ngradp1, bp);
         // std::cout << flux1;
         V qneg = dtpv*q;
         V qpos = dtpv*q;
@@ -312,7 +441,8 @@ int main()
             + ADB::constant(dtpv, bp)*(ops.div*flux1)
             - qtr_ad;
         res_norm = transport_residual.value().matrix().norm();
-        std::cout << "res_norm = " << res_norm << std::endl;
+        std::cout << "res_norm[" << it << "] = "
+                  << res_norm << std::endl;
 
         matr = transport_residual.derivative()[0];
         matr.makeCompressed();
@@ -329,10 +459,14 @@ int main()
         }
         // std::cout << x << std::endl;
         s1 = s.value() - x.array();
-        std::cerr << "Solve for s " << clock.secsSinceLast() << std::endl;
+        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]));
         }
-        std::cout << "s1 = \n" << s1 << std::endl;
+
+        it += 1;
     } while (res_norm > 1e-7);
+
+    std::cout << "Saturation solution:\ns1 = [\n" << s1 << "\n]\n";
 }