From bd240a31cb10aab7731d2f652f0d78cc560b6a31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 3 May 2013 16:24:20 +0200 Subject: [PATCH 1/9] Declare operator?() as 'const'. None of operator+(), operator-(), operator*() or operator/() modify their object. These methods are thus naturally declared as 'const'. --- AutoDiffBlock.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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); From d0866c39b5742869d16e5e3a14079ebb27aa58ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 3 May 2013 17:48:33 +0200 Subject: [PATCH 2/9] Implement first cut at upwind-mobility weighting for transport Specifically, - Introduce a new class, UpwindSelector which provides a selection service (method select()) to pick out a quantity across an interface in the upwind direction. Method select() is implemented only for the case of no counter-current flow (i.e., no effects of gravity and/or capillary pressure). The selector is affected through applying an explicit sparse matrix to each quantity. This may be inefficient, but only measurements will tell. - Split function fluxFunc() into two components, phaseMobility() and fluxFunc(). The former computes phase mobilities (per cell) for water and oil in a subset of the grid cells and the latter computes the flux function (for water, i.e., the first phase) based on a (2-) vector of phase mobilities--either per cell or per internal interface. - Compute explicit phase mobilities (per cell) in the non-linear loop before employing the upwind selector to pick out proper upwind mobilities per (internal) interface to define the quantity 'fw_face' which was previously computed by an arithmetic average. --- sim_simple.cpp | 104 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 91 insertions(+), 13 deletions(-) diff --git a/sim_simple.cpp b/sim_simple.cpp index 8e4128af8..cf87feffa 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,73 @@ struct HelperOps } }; -/// Returns fw(sw). +namespace { + 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 +207,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 +360,7 @@ int main() double res_norm = 1e100; V s1 = /*s0.leftCols<1>()*/0.5*V::Ones(nc,1); // Initial guess. + UpwindSelector upws(grid, ops); do { const std::vector& bp = block_pattern; ADB s = ADB::variable(0, s1, bp); @@ -296,9 +371,12 @@ int main() // 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::vector pmobc = phaseMobility(props, allcells, s.value()); + std::vector pmobf = upws.select(p1, pmobc); + + ADB fw_cell = fluxFunc(pmobc); + ADB fw_face = fluxFunc(pmobf); + // std::cout << fw_face; ADB flux1 = fw_face*ADB::constant(ngradp1, bp); // std::cout << flux1; From 9b65338461fe19b30a1c37e362f290451a44df2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 3 May 2013 20:21:37 +0200 Subject: [PATCH 3/9] Add MRST function for 'sim_simple' development --- mrst/sim_simple.m | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 mrst/sim_simple.m diff --git a/mrst/sim_simple.m b/mrst/sim_simple.m new file mode 100644 index 000000000..ef3943546 --- /dev/null +++ b/mrst/sim_simple.m @@ -0,0 +1,26 @@ +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]); + + state = initState(g, [], 0, repmat([ 0.5, 0.5 ], [g.cells.num, 1])); + state = incompTPFA(state, 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 From fc5a3f181e8933431b4a3fc0f4665c0e7b3d20ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 3 May 2013 20:28:36 +0200 Subject: [PATCH 4/9] Use S_o=1 as initial condition. This is the case in 'sim_simple.cpp', too. --- mrst/sim_simple.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/mrst/sim_simple.m b/mrst/sim_simple.m index ef3943546..12d615c08 100644 --- a/mrst/sim_simple.m +++ b/mrst/sim_simple.m @@ -14,8 +14,9 @@ function x = sim_simple(cartDims, physDims, tf, verb) src = addSource([], [1, g.cells.num], [1, -1], 'sat', [ 1, 0 ; 0, 1]); - state = initState(g, [], 0, repmat([ 0.5, 0.5 ], [g.cells.num, 1])); - state = incompTPFA(state, g, T, fluid, 'src', src); + 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, ... From 8b12d1506d65925f6a0ea773eefe251a1c3246bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 3 May 2013 23:45:13 +0200 Subject: [PATCH 5/9] Add MATLAB-compatible output of ADB::M values. This is mostly intended as an interim aid in development. Not enabled if 'NDEBUG' is defined. --- sim_simple.cpp | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/sim_simple.cpp b/sim_simple.cpp index cf87feffa..67fb8a378 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -112,7 +112,49 @@ struct HelperOps } }; +#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: From e31101e80fd1a3a7deef9fe1d3f10bed7fafeaea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 5 May 2013 16:44:51 +0200 Subject: [PATCH 6/9] Incorporate transmissibility into water flux definition. This means that we're computing the correct solution for dt=0.0005 (as verified by 'sim_simple.m', but we're still failing larger time steps). --- sim_simple.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/sim_simple.cpp b/sim_simple.cpp index 67fb8a378..965dfad4e 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -403,6 +403,8 @@ 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)); do { const std::vector& bp = block_pattern; ADB s = ADB::variable(0, s1, bp); @@ -411,16 +413,11 @@ 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; std::vector pmobc = phaseMobility(props, allcells, s.value()); std::vector pmobf = upws.select(p1, pmobc); ADB fw_cell = fluxFunc(pmobc); - ADB fw_face = fluxFunc(pmobf); - - // std::cout << fw_face; - ADB flux1 = fw_face*ADB::constant(ngradp1, bp); + ADB flux1 = pmobf[0] * nkdp; // std::cout << flux1; V qneg = dtpv*q; V qpos = dtpv*q; From 6f22d10b0678e778efb24a8aebd26d76134b65c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 5 May 2013 21:52:34 +0200 Subject: [PATCH 7/9] Do real splitting based on total flux/fractional flow We are now able to compute the correct solution with dt=0.5, but the convergence rate is dismal. --- sim_simple.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/sim_simple.cpp b/sim_simple.cpp index 965dfad4e..319add4e6 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -405,6 +405,11 @@ int main() 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 = (pmobf0[0] + pmobf0[1]) * ADB::function(nkdp.value(), null); do { const std::vector& bp = block_pattern; ADB s = ADB::variable(0, s1, bp); @@ -417,7 +422,8 @@ int main() std::vector pmobf = upws.select(p1, pmobc); ADB fw_cell = fluxFunc(pmobc); - ADB flux1 = pmobf[0] * nkdp; + const ADB fw_face = fluxFunc(pmobf); + ADB flux1 = fw_face * dflux; // std::cout << flux1; V qneg = dtpv*q; V qpos = dtpv*q; From 056b7ae29229780d36141b397b23684c885cebd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 5 May 2013 22:02:40 +0200 Subject: [PATCH 8/9] Make Darcy flux actually constant. The definition from commit 6f22d10 included a non-zero derivative which interfered with the automatic differentiation solver. Using a fully constant Darcy flux enables using dt=150 while still computing the correct solution. While here, output the solution (s1) using more precision/digits for comparison with 'sim_simple.m'. --- sim_simple.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/sim_simple.cpp b/sim_simple.cpp index 319add4e6..7373e164e 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -409,7 +409,12 @@ int main() 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 = (pmobf0[0] + pmobf0[1]) * ADB::function(nkdp.value(), null); + 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); + do { const std::vector& bp = block_pattern; ADB s = ADB::variable(0, s1, bp); From b06b0c99ca53ca31c2239e4fb59578e5c162a2f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 5 May 2013 22:57:42 +0200 Subject: [PATCH 9/9] Decrease amount of output and count iterations in the process. Specifically, output the solution only at the end (in MATLAB compatible format, for easy comparison) and annotate the diagnostic output (i.e., cpu time and residual norm) with the iteration number. --- sim_simple.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/sim_simple.cpp b/sim_simple.cpp index 7373e164e..f7ef084af 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -415,6 +415,7 @@ int main() 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); @@ -440,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(); @@ -457,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"; }