opm-flowdiagnostics updated to 80190c28ae0fd4a82cfe88c827559029982db83b

opm-flowdiagnostics-applications updated to 01ecb8fc22cb70d883edaad398bffc49878859c3
Which fixes error in TOF calculations.
This commit is contained in:
Jacob Støren 2017-01-05 14:57:25 +01:00
parent ebddaa31e6
commit 95206315fa
11 changed files with 391 additions and 42 deletions

View File

@ -13,5 +13,8 @@ set(ERT_GITHUB_SHA "236164870f011305aed2eca85c45944b021e4107")
set(OPM_COMMON_GITHUB_SHA "1216bc052542f24ec6fcfbe1947d52e6300ff754")
set(OPM_PARSER_GITHUB_SHA "a3496df501a4369fd827fbf0ff893d03deff425f")
set(OPM_FLOWDIAGNOSTICS_SHA "80190c28ae0fd4a82cfe88c827559029982db83b")
set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "01ecb8fc22cb70d883edaad398bffc49878859c3")
set(STRPRODUCTVER ${RESINSIGHT_MAJOR_VERSION}.${RESINSIGHT_MINOR_VERSION}.${RESINSIGHT_INCREMENT_VERSION})

View File

@ -30,7 +30,9 @@ list (APPEND TEST_SOURCE_FILES
)
list (APPEND EXAMPLE_SOURCE_FILES
examples/computeLocalSolutions.cpp
examples/computeToFandTracers.cpp
examples/computeTracers.cpp
)
list (APPEND PUBLIC_HEADER_FILES

View File

@ -0,0 +1,69 @@
/*
Copyright 2017 SINTEF ICT, Applied Mathematics.
Copyright 2017 Statoil ASA.
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/>.
*/
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include "exampleSetup.hpp"
#include <opm/flowdiagnostics/CellSet.hpp>
// Syntax (typical):
// computeToFandTracers case=<ecl_case_prefix> step=<report_number>
int main(int argc, char* argv[])
try {
example::Setup setup(argc, argv);
auto& fdTool = setup.toolbox;
// Create start sets from injector wells.
std::vector<Opm::FlowDiagnostics::CellSet> start;
for (const auto& well : setup.well_fluxes) {
if (!well.is_injector_well) {
continue;
}
std::vector<int> completion_cells;
completion_cells.reserve(well.completions.size());
for (const auto& completion : well.completions) {
const int grid_index = completion.grid_index;
const auto& ijk = completion.ijk;
const int cell_index = setup.graph.activeCell(ijk, grid_index);
if (cell_index >= 0) {
completion_cells.push_back(cell_index);
}
}
start.emplace_back(Opm::FlowDiagnostics::CellSetID(well.name), completion_cells);
}
// Solve for injection time of flight and tracers.
auto sol = fdTool.computeInjectionDiagnostics(start);
// Get local tof for first injector.
const auto& tof = sol.fd.timeOfFlight(start.front().id());
// Write it to standard out.
std::cout.precision(16);
for (auto item : tof) {
std::cout << item.first << " " << item.second << '\n';
}
}
catch (const std::exception& e) {
std::cerr << "Caught exception: " << e.what() << '\n';
}

View File

@ -29,7 +29,8 @@
// computeToFandTracers case=<ecl_case_prefix> step=<report_number>
int main(int argc, char* argv[])
try {
auto fdTool = example::setup(argc, argv);
example::Setup setup(argc, argv);
auto& fdTool = setup.toolbox;
// Solve for time of flight.
std::vector<Opm::FlowDiagnostics::CellSet> start;

View File

@ -0,0 +1,47 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
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/>.
*/
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include "exampleSetup.hpp"
int main(int argc, char** argv)
try {
example::Setup setup(argc, argv);
auto& fdTool = setup.toolbox;
// Solve for tracers.
Opm::FlowDiagnostics::CellSetID id("Example start set ID");
Opm::FlowDiagnostics::CellSet start(id, {123});
auto sol = fdTool.computeInjectionDiagnostics({start});
const auto& conc_sol = sol.fd.concentration(id);
// Write it to standard out.
std::cout.precision(16);
for (auto item : conc_sol) {
std::cout << item.first << " " << item.second << '\n';
}
}
catch (const std::exception& e) {
std::cerr << "Caught exception: " << e.what() << '\n';
}

View File

@ -158,32 +158,8 @@ namespace example {
}
}
inline Opm::FlowDiagnostics::Toolbox
initialiseFlowDiagnostics(const Opm::ECLGraph& G)
{
const auto connGraph = Opm::FlowDiagnostics::
ConnectivityGraph{ static_cast<int>(G.numCells()),
G.neighbours() };
// Create the Toolbox.
auto tool = Opm::FlowDiagnostics::Toolbox{ connGraph };
tool.assignPoreVolume(G.poreVolume());
tool.assignConnectionFlux(Hack::convert_flux_to_SI(extractFluxField(G)));
auto wsol = Opm::ECLWellSolution{};
const auto well_fluxes =
wsol.solution(G.rawResultData(), G.numGrids());
tool.assignInflowFlux(extractWellFlows(G, well_fluxes));
return tool;
}
inline auto setup(int argc, char* argv[])
-> decltype(initialiseFlowDiagnostics(std::declval<Opm::ECLGraph>()))
inline Opm::ECLGraph
initGraph(int argc, char* argv[])
{
// Obtain parameters from command line (possibly specifying a parameter file).
const bool verify_commandline_syntax = true;
@ -216,9 +192,51 @@ namespace example {
throw std::domain_error(os.str());
}
return initialiseFlowDiagnostics(graph);
return graph;
}
inline std::vector<Opm::ECLWellSolution::WellData>
initWellFluxes(const Opm::ECLGraph& G)
{
auto wsol = Opm::ECLWellSolution{};
return wsol.solution(G.rawResultData(), G.numGrids());
}
inline Opm::FlowDiagnostics::Toolbox
initToolbox(const Opm::ECLGraph& G, const std::vector<Opm::ECLWellSolution::WellData>& well_fluxes)
{
const auto connGraph = Opm::FlowDiagnostics::
ConnectivityGraph{ static_cast<int>(G.numCells()),
G.neighbours() };
// Create the Toolbox.
auto tool = Opm::FlowDiagnostics::Toolbox{ connGraph };
tool.assignPoreVolume(G.poreVolume());
tool.assignConnectionFlux(Hack::convert_flux_to_SI(extractFluxField(G)));
tool.assignInflowFlux(extractWellFlows(G, well_fluxes));
return tool;
}
struct Setup
{
Setup(int argc, char** argv)
: graph(initGraph(argc, argv))
, well_fluxes(initWellFluxes(graph))
, toolbox(initToolbox(graph, well_fluxes))
{
}
Opm::ECLGraph graph;
std::vector<Opm::ECLWellSolution::WellData> well_fluxes;
Opm::FlowDiagnostics::Toolbox toolbox;
};
} // namespace example

View File

@ -202,8 +202,9 @@ namespace Opm
// Otherwise: add data for this well.
WellData wd;
wd.name = trimSpacesRight(zwel[well * ih.nzwel]);
const int ncon = iwel[well * ih.niwel + IWEL_CONNECTIONS_INDEX];
const bool is_producer = (iwel[well * ih.niwel + IWEL_TYPE_INDEX] == IWEL_TYPE_PRODUCER);
wd.is_injector_well = !is_producer;
const int ncon = iwel[well * ih.niwel + IWEL_CONNECTIONS_INDEX];
wd.completions.reserve(ncon);
for (int comp_index = 0; comp_index < ncon; ++comp_index) {
const int icon_offset = (well*ih.ncwma + comp_index) * ih.nicon;

View File

@ -139,7 +139,7 @@ Toolbox::Impl::injDiag(const std::vector<CellSet>& start_sets)
using Conc = Solution::TracerConcentration;
TracerTofSolver solver(downstream_conn_, upstream_conn_, pvol_, only_inflow_flux_);
sol.assignGlobalToF(solver.solveGlobal(start_sets));
sol.assignGlobalToF(solver.solveGlobal());
for (const auto& start : start_sets) {
auto solution = solver.solveLocal(start);
@ -167,7 +167,7 @@ Toolbox::Impl::prodDiag(const std::vector<CellSet>& start_sets)
using Conc = Solution::TracerConcentration;
TracerTofSolver solver(upstream_conn_, downstream_conn_, pvol_, only_outflow_flux_);
sol.assignGlobalToF(solver.solveGlobal(start_sets));
sol.assignGlobalToF(solver.solveGlobal());
for (const auto& start : start_sets) {
auto solution = solver.solveLocal(start);
@ -198,7 +198,7 @@ Toolbox::Impl::buildAssembledConnections()
if (connection_flux > 0.0) {
downstream_conn_.addConnection(cells.first, cells.second, connection_flux);
upstream_conn_.addConnection(cells.second, cells.first, connection_flux);
} else {
} else if (connection_flux < 0.0) {
downstream_conn_.addConnection(cells.second, cells.first, -connection_flux);
upstream_conn_.addConnection(cells.first, cells.second, -connection_flux);
}

View File

@ -107,13 +107,11 @@ namespace FlowDiagnostics
std::vector<double> TracerTofSolver::solveGlobal(const std::vector<CellSet>& all_startsets)
std::vector<double> TracerTofSolver::solveGlobal()
{
// Reset solver variables and set source terms.
prepareForSolve();
for (const CellSet& startset : all_startsets) {
setupStartArray(startset);
}
setupStartArrayFromSource();
// Compute topological ordering and solve.
computeOrdering();
@ -183,6 +181,20 @@ namespace FlowDiagnostics
void TracerTofSolver::setupStartArrayFromSource()
{
const int num_cells = pv_.size();
for (int cell = 0; cell < num_cells; ++cell) {
if (source_term_[cell] > 0.0) {
is_start_[cell] = 1;
}
}
}
void TracerTofSolver::computeOrdering()
{
// Compute reverse topological ordering.
@ -239,7 +251,7 @@ namespace FlowDiagnostics
}
// Extract data from solution.
sequence_.resize(num_cells);
sequence_.resize(num_cells); // For local solutions this is the upper limit of the size. TODO: use exact size.
const int num_comp = tarjan_get_numcomponents(result.get());
component_starts_.resize(num_comp + 1);
component_starts_[0] = 0;
@ -305,13 +317,13 @@ namespace FlowDiagnostics
for (const auto& conn : g_reverse_.cellNeighbourhood(cell)) {
const int upwind_cell = conn.neighbour;
const double flux = conn.weight;
upwind_tof_contrib += tof_[upwind_cell] * flux;
upwind_tof_contrib += tof_[upwind_cell] * tracer_[upwind_cell] * flux;
upwind_tracer_contrib += tracer_[upwind_cell] * flux;
}
// Compute time-of-flight and tracer.
tof_[cell] = (pv_[cell] + upwind_tof_contrib) / total_influx;
tracer_[cell] = is_start_[cell] ? 1.0 : upwind_tracer_contrib / total_influx;
tof_[cell] = (pv_[cell]*tracer_[cell] + upwind_tof_contrib) / (total_influx*tracer_[cell]);
}

View File

@ -59,7 +59,7 @@ namespace FlowDiagnostics
/// Compute the global (combining all sources) time-of-flight of each cell.
///
/// TODO: also compute tracer solution.
std::vector<double> solveGlobal(const std::vector<CellSet>& all_startsets);
std::vector<double> solveGlobal();
/// Output data struct for solveLocal().
struct LocalSolution {
@ -111,6 +111,8 @@ namespace FlowDiagnostics
void setupStartArray(const CellSet& startset);
void setupStartArrayFromSource();
void computeOrdering();
void computeLocalOrdering(const CellSet& startset);

View File

@ -298,7 +298,7 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
const auto tof = fwd.fd.timeOfFlight();
BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
std::vector<double> expected = { 1.0, 2.0, 3.0, 4.0, 0.0 };
std::vector<double> expected = { 1.0, 2.0, 3.0, 4.0, 5.0 };
check_is_close(tof, expected);
}
@ -307,7 +307,7 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
const auto tof = rev.fd.timeOfFlight();
BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
std::vector<double> expected = { 0.0, 4.0, 3.0, 2.0, 1.0 };
std::vector<double> expected = { 5.0, 4.0, 3.0, 2.0, 1.0 };
check_is_close(tof, expected);
}
@ -439,4 +439,198 @@ BOOST_AUTO_TEST_CASE (OneDimCase)
}
// Arrows indicate a flux of 0.3, O is a source of 0.3
// and X is a sink of 0.3 (each cell has a pore volume of 0.3).
// ----------------------------
// | | | |
// | O -> -> |
// | | -> |
// | | | || |
// -------------^--------VV----
// | | | | |
// | | | |
// | O -> | XX |
// | | | |
// ----------------------------
// Cell indices:
// ----------------------------
// | | | |
// | | | |
// | 3 | 4 | 5 |
// | | | |
// ----------------------------
// | | | |
// | | | |
// | 0 | 1 | 2 |
// | | | |
// ----------------------------
// Expected global injection TOF:
// ----------------------------
// | | | |
// | | | |
// | 1.0 | 2.0 | 2.5 |
// | | | |
// ----------------------------
// | | | |
// | | | |
// | 1.0 | 2.0 | 3.0 |
// | | | |
// ----------------------------
// Expected global production TOF:
// ----------------------------
// | | | |
// | | | |
// | 2.5 | 1.5 | 1.0 |
// | | | |
// ----------------------------
// | | | |
// | | | |
// | 3.5 | 2.5 | 0.5 |
// | | | |
// ----------------------------
BOOST_AUTO_TEST_CASE (LocalSolutions)
{
using namespace Opm::FlowDiagnostics;
const auto cas = Setup(3, 2);
const auto& graph = cas.connectivity();
// Create fluxes.
ConnectionValues flux(ConnectionValues::NumConnections{ graph.numConnections() },
ConnectionValues::NumPhases { 1 });
const size_t nconn = cas.connectivity().numConnections();
for (size_t conn = 0; conn < nconn; ++conn) {
BOOST_TEST_MESSAGE("Connection " << conn << " connects cells "
<< graph.connection(conn).first << " and "
<< graph.connection(conn).second);
}
using C = ConnectionValues::ConnID;
using P = ConnectionValues::PhaseID;
flux(C{0}, P{0}) = 0.3;
flux(C{1}, P{0}) = 0.0;
flux(C{2}, P{0}) = 0.3;
flux(C{3}, P{0}) = 0.6;
flux(C{4}, P{0}) = 0.0;
flux(C{5}, P{0}) = 0.3;
flux(C{6}, P{0}) = -0.6;
// Create well in/out flows.
CellSetValues wellflow = { {0, 0.3}, {3, 0.3}, {2, -0.6} };
Toolbox diagTool(graph);
diagTool.assignPoreVolume(cas.poreVolume());
diagTool.assignConnectionFlux(flux);
diagTool.assignInflowFlux(wellflow);
auto injstart = std::vector<CellSet>{ CellSet(CellSetID("I-1"), {0}),
CellSet(CellSetID("I-2"), {3}) };
auto prdstart = std::vector<CellSet>{ CellSet(CellSetID("P-1"), {2}) };
const auto fwd = diagTool.computeInjectionDiagnostics(injstart);
const auto rev = diagTool.computeProductionDiagnostics(prdstart);
// Global ToF field (accumulated from all injectors)
{
const auto tof = fwd.fd.timeOfFlight();
BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
std::vector<double> expected = { 1.0, 2.0, 3.0, 1.0, 2.0, 2.5 };
check_is_close(tof, expected);
}
// Global ToF field (accumulated from all producers)
{
const auto tof = rev.fd.timeOfFlight();
BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
std::vector<double> expected = { 3.5, 2.5, 0.5, 2.5, 1.5, 1.0 };
check_is_close(tof, expected);
}
// Verify set of start points.
{
using VCS = std::vector<Opm::FlowDiagnostics::CellSet>;
using VCSI = std::vector<Opm::FlowDiagnostics::CellSetID>;
using P = std::pair<VCS, VCSI>;
std::vector<P> pairs { P{ injstart, fwd.fd.startPoints() }, P{ prdstart, rev.fd.startPoints() } };
for (const auto& p : pairs) {
const auto& s1 = p.first;
const auto& s2 = p.second;
BOOST_CHECK_EQUAL(s1.size(), s2.size());
for (const auto& pt : s2) {
// ID of 'pt' *MUST* be in set of identified start points.
auto pos = std::find_if(s1.begin(), s1.end(),
[&pt](const CellSet& s)
{
return s.id().to_string() == pt.to_string();
});
BOOST_CHECK(pos != s1.end());
}
}
}
// Local I-1 tracer concentration.
{
const auto conc = fwd.fd.concentration(CellSetID("I-1"));
std::vector<std::pair<int, double>> expected = { {0, 1.0}, {1, 1.0}, {2, 0.5}, {4, 0.5}, {5, 0.5} };
BOOST_REQUIRE_EQUAL(conc.size(), expected.size());
int i = 0;
for (const auto& v : conc) {
BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
BOOST_CHECK_EQUAL(v.first, expected[i].first);
BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
++i;
}
}
// Local I-1 tof.
{
const auto tof = fwd.fd.timeOfFlight(CellSetID("I-1"));
std::vector<std::pair<int, double>> expected = { {0, 1.0}, {1, 2.0}, {2, 3.5}, {4, 2.5}, {5, 3.0} };
BOOST_REQUIRE_EQUAL(tof.size(), expected.size());
int i = 0;
for (const auto& v : tof) {
BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
BOOST_CHECK_EQUAL(v.first, expected[i].first);
BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
++i;
}
}
// Local I-2 tracer concentration.
{
const auto conc = fwd.fd.concentration(CellSetID("I-2"));
std::vector<std::pair<int, double>> expected = { {2, 0.5}, {3, 1.0}, {4, 0.5}, {5, 0.5} };
BOOST_REQUIRE_EQUAL(conc.size(), expected.size());
int i = 0;
for (const auto& v : conc) {
BOOST_TEST_MESSAGE("Conc[" << v.first << "] = " << v.second);
BOOST_CHECK_EQUAL(v.first, expected[i].first);
BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
++i;
}
}
// Local I-2 tof.
{
const auto tof = fwd.fd.timeOfFlight(CellSetID("I-2"));
std::vector<std::pair<int, double>> expected = { {2, 2.5}, {3, 1.0}, {4, 1.5}, {5, 2.0} };
BOOST_REQUIRE_EQUAL(tof.size(), expected.size());
int i = 0;
for (const auto& v : tof) {
BOOST_TEST_MESSAGE("ToF[" << v.first << "] = " << v.second);
BOOST_CHECK_EQUAL(v.first, expected[i].first);
BOOST_CHECK_CLOSE(v.second, expected[i].second, 1.0e-10);
++i;
}
}
}
BOOST_AUTO_TEST_SUITE_END()