OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added calculation of flow diagnostics quantities.

Browse Source 
New functions:
 - computeFandPhi()
 - computeLorenz()
 - computeSweep()

Also a unit test has been added for the new features.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-01-21 14:58:44 +01:00
parent 0f6d2104d4
commit 1e0d2ec43e
4 changed files with 448 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -112,6 +112,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/simulator/SimulatorTimer.cpp
opm/core/tof/AnisotropicEikonal.cpp
opm/core/tof/DGBasis.cpp
opm/core/tof/FlowDiagnostics.cpp
opm/core/tof/TofReorder.cpp
opm/core/tof/TofDiscGalReorder.cpp
opm/core/transport/TransportSolverTwophaseInterface.cpp
@ -164,6 +165,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_ug.cpp
tests/test_cubic.cpp
tests/test_event.cpp
tests/test_flowdiagnostics.cpp
tests/test_nonuniformtablelinear.cpp
tests/test_sparsevector.cpp
tests/test_sparsetable.cpp
@ -383,6 +385,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/initStateEquil_impl.hpp
opm/core/tof/AnisotropicEikonal.hpp
opm/core/tof/DGBasis.hpp
opm/core/tof/FlowDiagnostics.hpp
opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp
opm/core/transport/TransportSolverTwophaseInterface.hpp

165
opm/core/tof/FlowDiagnostics.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/tof/FlowDiagnostics.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <numeric>
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<double>, std::vector<double>> computeFandPhi(const std::vector<double>& pv,
const std::vector<double>& ftof,
const std::vector<double>& 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<double, double> D2;
std::vector<D2> 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<double> 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<double> 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<double>& flowcap,
const std::vector<double>& 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<double>, std::vector<double>> computeSweep(const std::vector<double>& flowcap,
const std::vector<double>& 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<double> Ev;
std::vector<double> 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

82
opm/core/tof/FlowDiagnostics.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED
#define OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED
#include <vector>
#include <utility>
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<double>, std::vector<double>> computeFandPhi(const std::vector<double>& pv,
const std::vector<double>& ftof,
const std::vector<double>& 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<double>& flowcap,
const std::vector<double>& 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<double>, std::vector<double>> computeSweep(const std::vector<double>& flowcap,
const std::vector<double>& storagecap);
} // namespace Opm
#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED

198
tests/test_flowdiagnostics.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <opm/core/tof/FlowDiagnostics.hpp>
const std::vector<double> pv(16, 18750.0);
const std::vector<double> 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<double> 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<double> 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<double> 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<double> 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<double> 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<double> wrong_length(7, 0.0);
using namespace Opm;
template <class C>
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);
}