diff --git a/AutoDiffBlock.hpp b/AutoDiffBlock.hpp index ff5041347..4a8fe2706 100644 --- a/AutoDiffBlock.hpp +++ b/AutoDiffBlock.hpp @@ -76,7 +76,7 @@ namespace AutoDiff } /// Operator + - ForwardBlock operator+(const ForwardBlock& rhs) + ForwardBlock operator+(const ForwardBlock& rhs) const { std::vector 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 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 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 jac(num_blocks); diff --git a/mrst/sim_simple.m b/mrst/sim_simple.m new file mode 100644 index 000000000..12d615c08 --- /dev/null +++ b/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 diff --git a/sim_simple.cpp b/sim_simple.cpp index 8e4128af8..f7ef084af 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -55,7 +55,8 @@ struct HelperOps typedef AutoDiff::ForwardBlock::V V; /// A list of internal faces. - Eigen::Array internal_faces; + typedef Eigen::Array 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 +#endif // !defined(NDEBUG) + +namespace { +#if !defined(NDEBUG) + void + printSparseMatrix(const Eigen::SparseMatrix& A, + std::FILE* fp) + { + typedef Eigen::SparseMatrix::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(i[k] + 1), + static_cast(j + 1), x[k]); + } + } + } + + void + printSparseMatrix(const Eigen::SparseMatrix& 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 + class UpwindSelector { + public: + typedef AutoDiff::ForwardBlock 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 + select(const typename ADB::V& press, + const std::vector& xc ) const + { + typedef HelperOps::iFaces::Index ifIndex; + const ifIndex nif = h_.internal_faces.size(); + + // Define selector structure. + typedef typename Eigen::Triplet Triplet; + std::vector 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 xf; xf.reserve(xc.size()); + for (typename std::vector::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 -ADB fluxFunc(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const typename ADB::V& sw) +std::vector +phaseMobility(const Opm::IncompPropertiesInterface& props, + const std::vector& cells, + const typename ADB::V& sw) { typedef Eigen::Array TwoCol; typedef Eigen::Array FourCol; @@ -144,11 +249,22 @@ ADB fluxFunc(const Opm::IncompPropertiesInterface& props, const double* mu = props.viscosity(); std::vector dmw = { krwjac/mu[0] }; std::vector dmo = { krojac/mu[1] }; - ADB mw_ad = ADB::function(krw/mu[0], dmw); - ADB mo_ad = ADB::function(kro/mu[1], dmo); - ADB fw = mw_ad / (mw_ad + mo_ad); - // std::cout << mw_ad << mo_ad << (mw_ad + mo_ad) << fw; - return fw; + + std::vector pmobc = { ADB::function(krw / mu[0], dmw) , + ADB::function(kro / mu[1], dmo) }; + return pmobc; +} + +/// Returns fw(sw). +template +ADB +fluxFunc(const std::vector& 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 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 pmobc0 = phaseMobility(props, allcells, s00.value()); + const std::vector pmobf0 = upws.select(p1, pmobc0); + const std::vector 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& bp = block_pattern; ADB s = ADB::variable(0, s1, bp); @@ -294,13 +424,12 @@ int main() * Eigen::Map(grid.cell_volumes, nc, 1); V dtpv = dt/pv; // std::cout << dtpv; - V ngradp1 = ops.ngrad*p1.matrix(); - // std::cout << ngradp1 << std::endl; - ADB fw_cell = fluxFunc(props, allcells, s.value()); - // std::cout << fw_cell; - ADB fw_face = ops.caver*fw_cell; - // std::cout << fw_face; - ADB flux1 = fw_face*ADB::constant(ngradp1, bp); + std::vector pmobc = phaseMobility(props, allcells, s.value()); + std::vector pmobf = upws.select(p1, pmobc); + + 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; @@ -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"; }