From 77ae151be909d7a43b22bee38e29ef172420b8e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 21 Jan 2015 14:58:44 +0100 Subject: [PATCH] Added calculation of flow diagnostics quantities. New functions: - computeFandPhi() - computeLorenz() - computeSweep() Also a unit test has been added for the new features. --- opm/core/tof/FlowDiagnostics.cpp | 165 ++++++++++++++++++++++++++ opm/core/tof/FlowDiagnostics.hpp | 82 +++++++++++++ tests/test_flowdiagnostics.cpp | 198 +++++++++++++++++++++++++++++++ 3 files changed, 445 insertions(+) create mode 100644 opm/core/tof/FlowDiagnostics.cpp create mode 100644 opm/core/tof/FlowDiagnostics.hpp create mode 100644 tests/test_flowdiagnostics.cpp diff --git a/opm/core/tof/FlowDiagnostics.cpp b/opm/core/tof/FlowDiagnostics.cpp new file mode 100644 index 000000000..ef3a4bc6d --- /dev/null +++ b/opm/core/tof/FlowDiagnostics.cpp @@ -0,0 +1,165 @@ +/* + Copyright 2015 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 . +*/ + +#include + +#include +#include + +namespace Opm +{ + + + /// \brief Compute flow-capacity/storage-capacity based on time-of-flight. + /// + /// The F-Phi curve is an analogue to the fractional flow curve in a 1D + /// displacement. It can be used to compute other interesting diagnostic + /// quantities such as the Lorenz coefficient. For a technical description + /// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625). + /// + /// \param[in] pv pore volumes of each cell + /// \param[in] ftof forward (time from injector) time-of-flight values for each cell + /// \param[in] rtof reverse (time to producer) time-of-flight values for each cell + /// \return a pair of vectors, the first containing F (flow capacity) the second + /// containing Phi (storage capacity). + std::pair, std::vector> computeFandPhi(const std::vector& pv, + const std::vector& ftof, + const std::vector& rtof) + { + if (pv.size() != ftof.size() || pv.size() != rtof.size()) { + OPM_THROW(std::runtime_error, "computeFandPhi(): Input vectors must have same size."); + } + + // Sort according to total travel time. + const int n = pv.size(); + typedef std::pair D2; + std::vector time_and_pv(n); + for (int ii = 0; ii < n; ++ii) { + time_and_pv[ii].first = ftof[ii] + rtof[ii]; // Total travel time. + time_and_pv[ii].second = pv[ii]; + } + std::sort(time_and_pv.begin(), time_and_pv.end()); + + // Compute Phi. + std::vector Phi(n + 1); + Phi[0] = 0.0; + for (int ii = 0; ii < n; ++ii) { + Phi[ii+1] = time_and_pv[ii].second; + } + std::partial_sum(Phi.begin(), Phi.end(), Phi.begin()); + const double vt = Phi.back(); // Total pore volume. + for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop. + Phi[ii] /= vt; // Normalize Phi. + } + + // Compute F. + std::vector F(n + 1); + F[0] = 0.0; + for (int ii = 0; ii < n; ++ii) { + F[ii+1] = time_and_pv[ii].second / time_and_pv[ii].first; + } + std::partial_sum(F.begin(), F.end(), F.begin()); + const double ft = F.back(); // Total flux. + for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop. + F[ii] /= ft; // Normalize Phi. + } + + return std::make_pair(F, Phi); + } + + + + + + /// \brief Compute the Lorenz coefficient based on the F-Phi curve. + /// + /// The Lorenz coefficient is a measure of heterogeneity. It is equal + /// to twice the area between the F-Phi curve and the F = Phi line. + /// The coefficient can vary from zero to one. If the coefficient is + /// zero (so the F-Phi curve is a straight line) we have perfect + /// piston-like displacement while a coefficient of one indicates + /// infinitely heterogenous displacement (essentially no sweep). + /// + /// Note: The coefficient is analogous to the Gini coefficient of + /// economic theory, where the name Lorenz curve is applied to + /// what we call the F-Phi curve. + /// + /// \param[in] flowcap flow capacity (F) as from computeFandPhi() + /// \param[in] storagecap storage capacity (Phi) as from computeFandPhi() + /// \return the Lorenz coefficient + double computeLorenz(const std::vector& flowcap, + const std::vector& storagecap) + { + if (flowcap.size() != storagecap.size()) { + OPM_THROW(std::runtime_error, "computeLorenz(): Input vectors must have same size."); + } + double integral = 0.0; + // Trapezoid quadrature of the curve F(Phi). + const int num_intervals = flowcap.size() - 1; + for (int ii = 0; ii < num_intervals; ++ii) { + const double len = storagecap[ii+1] - storagecap[ii]; + integral += (flowcap[ii] + flowcap[ii+1]) * len / 2.0; + } + return 2.0 * (integral - 0.5); + } + + + + + + /// \brief Compute sweep efficiency versus dimensionless time (PVI). + /// + /// The sweep efficiency is analogue to 1D displacement using the + /// F-Phi curve as flux function. + /// + /// \param[in] flowcap flow capacity (F) as from computeFandPhi() + /// \param[in] storagecap storage capacity (Phi) as from computeFandPhi() + /// \return a pair of vectors, the first containing Ev (sweep efficiency) + /// the second containing tD (dimensionless time). + std::pair, std::vector> computeSweep(const std::vector& flowcap, + const std::vector& storagecap) + { + if (flowcap.size() != storagecap.size()) { + OPM_THROW(std::runtime_error, "computeSweep(): Input vectors must have same size."); + } + + // Compute tD and Ev simultaneously, + // skipping identical Phi data points. + const int n = flowcap.size(); + std::vector Ev; + std::vector tD; + tD.reserve(n); + Ev.reserve(n); + tD.push_back(0.0); + Ev.push_back(0.0); + for (int ii = 1; ii < n; ++ii) { // Note loop limits. + const double fd = flowcap[ii] - flowcap[ii-1]; + const double sd = storagecap[ii] - storagecap[ii-1]; + if (fd != 0.0) { + tD.push_back(sd/fd); + Ev.push_back(storagecap[ii] + (1.0 - flowcap[ii]) * tD.back()); + } + } + + return std::make_pair(Ev, tD); + } + + + +} // namespace Opm diff --git a/opm/core/tof/FlowDiagnostics.hpp b/opm/core/tof/FlowDiagnostics.hpp new file mode 100644 index 000000000..faaff1ee5 --- /dev/null +++ b/opm/core/tof/FlowDiagnostics.hpp @@ -0,0 +1,82 @@ +/* + Copyright 2015 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 . +*/ + +#ifndef OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED +#define OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED + + +#include +#include + + +namespace Opm +{ + + /// \brief Compute flow-capacity/storage-capacity based on time-of-flight. + /// + /// The F-Phi curve is an analogue to the fractional flow curve in a 1D + /// displacement. It can be used to compute other interesting diagnostic + /// quantities such as the Lorenz coefficient. For a technical description + /// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625). + /// + /// \param[in] pv pore volumes of each cell + /// \param[in] ftof forward (time from injector) time-of-flight values for each cell + /// \param[in] rtof reverse (time to producer) time-of-flight values for each cell + /// \return a pair of vectors, the first containing F (flow capacity) the second + /// containing Phi (storage capacity). + std::pair, std::vector> computeFandPhi(const std::vector& pv, + const std::vector& ftof, + const std::vector& rtof); + + + /// \brief Compute the Lorenz coefficient based on the F-Phi curve. + /// + /// The Lorenz coefficient is a measure of heterogeneity. It is equal + /// to twice the area between the F-Phi curve and the F = Phi line. + /// The coefficient can vary from zero to one. If the coefficient is + /// zero (so the F-Phi curve is a straight line) we have perfect + /// piston-like displacement while a coefficient of one indicates + /// infinitely heterogenous displacement (essentially no sweep). + /// + /// Note: The coefficient is analogous to the Gini coefficient of + /// economic theory, where the name Lorenz curve is applied to + /// what we call the F-Phi curve. + /// + /// \param[in] flowcap flow capacity (F) as from computeFandPhi() + /// \param[in] storagecap storage capacity (Phi) as from computeFandPhi() + /// \return the Lorenz coefficient + double computeLorenz(const std::vector& flowcap, + const std::vector& storagecap); + + + /// \brief Compute sweep efficiency versus dimensionless time (PVI). + /// + /// The sweep efficiency is analogue to 1D displacement using the + /// F-Phi curve as flux function. + /// + /// \param[in] flowcap flow capacity (F) as from computeFandPhi() + /// \param[in] storagecap storage capacity (Phi) as from computeFandPhi() + /// \return a pair of vectors, the first containing Ev (sweep efficiency) + /// the second containing tD (dimensionless time). + std::pair, std::vector> computeSweep(const std::vector& flowcap, + const std::vector& storagecap); + +} // namespace Opm + +#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED diff --git a/tests/test_flowdiagnostics.cpp b/tests/test_flowdiagnostics.cpp new file mode 100644 index 000000000..a3d700e32 --- /dev/null +++ b/tests/test_flowdiagnostics.cpp @@ -0,0 +1,198 @@ +/* + Copyright 2015 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 . +*/ + +#include + +#if defined(HAVE_DYNAMIC_BOOST_TEST) +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + + +#define BOOST_TEST_MODULE FlowDiagnosticsTests +#include +#include + +const std::vector pv(16, 18750.0); + +const std::vector ftof = { + 5.399999999999992e+04, + 1.139999999999997e+05, + 2.819999999999993e+05, + 8.220000000000012e+05, + 1.139999999999998e+05, + 1.774285714285711e+05, + 3.160150375939849e+05, + 8.156820645778908e+05, + 2.819999999999994e+05, + 3.160150375939841e+05, + 3.935500938781204e+05, + 7.765612369042073e+05, + 8.220000000000000e+05, + 8.156820645778894e+05, + 7.765612369042063e+05, + 8.218220225906991e+05 +}; + +const std::vector rtof = { + 8.218220225906990e+05, + 7.765612369042046e+05, + 8.156820645778881e+05, + 8.219999999999976e+05, + 7.765612369042051e+05, + 3.935500938781204e+05, + 3.160150375939846e+05, + 2.820000000000001e+05, + 8.156820645778885e+05, + 3.160150375939850e+05, + 1.774285714285714e+05, + 1.140000000000000e+05, + 8.219999999999980e+05, + 2.819999999999998e+05, + 1.140000000000000e+05, + 5.400000000000003e+04 +}; + +const std::vector F = { + 0, + 9.568875799840706e-02, + 1.913775159968141e-01, + 2.778231480508526e-01, + 3.642687801048911e-01, + 4.266515906731506e-01, + 4.890344012414101e-01, + 5.503847464649610e-01, + 6.117350916885119e-01, + 6.730854369120627e-01, + 7.344357821356134e-01, + 7.842099754925904e-01, + 8.339841688495674e-01, + 8.837583622065442e-01, + 9.335325555635212e-01, + 9.667662777817606e-01, + 1.000000000000000e+00 +}; + +const std::vector Phi = { + 0, + 6.250000000000000e-02, + 1.250000000000000e-01, + 1.875000000000000e-01, + 2.500000000000000e-01, + 3.125000000000000e-01, + 3.750000000000000e-01, + 4.375000000000000e-01, + 5.000000000000000e-01, + 5.625000000000000e-01, + 6.250000000000000e-01, + 6.875000000000000e-01, + 7.500000000000000e-01, + 8.125000000000000e-01, + 8.750000000000000e-01, + 9.375000000000000e-01, + 1.000000000000000e+00 +}; + +const std::vector Ev = { + 0, + 6.531592770912591e-01, + 6.531592770912593e-01, + 7.096322601771997e-01, + 7.096322601772002e-01, + 8.869254748464411e-01, + 8.869254748464422e-01, + 8.955406718746977e-01, + 8.955406718746983e-01, + 8.955406718746991e-01, + 8.955406718746991e-01, + 9.584612275378565e-01, + 9.584612275378565e-01, + 9.584612275378569e-01, + 9.584612275378566e-01, + 1.000000000000000e+00, + 1.000000000000000e+00 +}; + +const std::vector tD = { + 0, + 6.531592770912591e-01, + 6.531592770912593e-01, + 7.229977792392133e-01, + 7.229977792392139e-01, + 1.001878553253259e+00, + 1.001878553253261e+00, + 1.018739173712224e+00, + 1.018739173712226e+00, + 1.018739173712227e+00, + 1.018739173712227e+00, + 1.255670776053656e+00, + 1.255670776053656e+00, + 1.255670776053659e+00, + 1.255670776053656e+00, + 1.880619919417231e+00, + 1.880619919417231e+00 +}; + +std::vector wrong_length(7, 0.0); + + + +using namespace Opm; + + +template +void compareCollections(const C& c1, const C& c2, const double tolerance = 1e-11) +{ + BOOST_REQUIRE(c1.size() == c2.size()); + auto c1it = c1.begin(); + auto c2it = c2.begin(); + for (; c1it != c1.end(); ++c1it, ++c2it) { + BOOST_CHECK_CLOSE(*c1it, *c2it, tolerance); + } +} + + +BOOST_AUTO_TEST_CASE(FandPhi) +{ + BOOST_CHECK_THROW(computeFandPhi(pv, ftof, wrong_length), std::runtime_error); + auto FPhi = computeFandPhi(pv, ftof, rtof); + compareCollections(FPhi.first, F); + compareCollections(FPhi.second, Phi); +} + + + + +BOOST_AUTO_TEST_CASE(Lorenz) +{ + BOOST_CHECK_THROW(computeLorenz(F, wrong_length), std::runtime_error); + const double Lc = computeLorenz(F, Phi); + BOOST_CHECK_CLOSE(Lc, 1.645920738950826e-01, 1e-11); +} + + + + +BOOST_AUTO_TEST_CASE(Sweep) +{ + BOOST_CHECK_THROW(computeSweep(F, wrong_length), std::runtime_error); + auto et = computeSweep(F, Phi); + compareCollections(et.first, Ev); + compareCollections(et.second, tD); +}