Updated unit test of flow-diag to allign with the example in the flowdiag-app library.

This commit is contained in:
Jacob Støren 2016-12-08 17:09:03 +01:00
parent 147cb5ebe0
commit e25f6cb84d
4 changed files with 178 additions and 188 deletions

View File

@ -0,0 +1,154 @@
/*
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/>.
*/
#ifndef OPM_EXAMPLESETUP_HEADER_INCLUDED
#define OPM_EXAMPLESETUP_HEADER_INCLUDED
#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
#include <opm/flowdiagnostics/ConnectionValues.hpp>
#include <opm/flowdiagnostics/Toolbox.hpp>
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLWellSolution.hpp>
#include <exception>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#include <boost/filesystem.hpp>
namespace example {
inline Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G)
{
using ConnVals = Opm::FlowDiagnostics::ConnectionValues;
using NConn = ConnVals::NumConnections;
using NPhas = ConnVals::NumPhases;
const auto nconn = NConn{G.numConnections()};
const auto nphas = NPhas{3};
auto flux = ConnVals(nconn, nphas);
auto phas = ConnVals::PhaseID{0};
for (const auto& p : { Opm::ECLGraph::PhaseIndex::Aqua ,
Opm::ECLGraph::PhaseIndex::Liquid ,
Opm::ECLGraph::PhaseIndex::Vapour })
{
const auto pflux = G.flux(p);
if (! pflux.empty()) {
assert (pflux.size() == nconn.total);
auto conn = ConnVals::ConnID{0};
for (const auto& v : pflux) {
flux(conn, phas) = v;
conn.id += 1;
}
}
phas.id += 1;
}
return flux;
}
template <class WellFluxes>
Opm::FlowDiagnostics::CellSetValues
extractWellFlows(const Opm::ECLGraph& G,
const WellFluxes& well_fluxes)
{
Opm::FlowDiagnostics::CellSetValues inflow;
for (const auto& well : well_fluxes) {
for (const auto& completion : well.completions) {
const int grid_index = completion.grid_index;
const auto& ijk = completion.ijk;
const int cell_index = G.activeCell(ijk, grid_index);
if (cell_index >= 0) {
inflow.emplace(cell_index, completion.reservoir_inflow_rate);
}
}
}
return inflow;
}
namespace Hack {
inline Opm::FlowDiagnostics::ConnectionValues
convert_flux_to_SI(Opm::FlowDiagnostics::ConnectionValues&& fl)
{
using Co = Opm::FlowDiagnostics::ConnectionValues::ConnID;
using Ph = Opm::FlowDiagnostics::ConnectionValues::PhaseID;
const auto nconn = fl.numConnections();
const auto nphas = fl.numPhases();
for (auto phas = Ph{0}; phas.id < nphas; ++phas.id) {
for (auto conn = Co{0}; conn.id < nconn; ++conn.id) {
fl(conn, phas) /= 86400;
}
}
return fl;
}
}
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;
}
} // namespace example
#endif // OPM_EXAMPLESETUP_HEADER_INCLUDED

View File

@ -1,131 +1,45 @@
#include "gtest/gtest.h"
#include <opm/flowdiagnostics/CellSet.hpp>
#include <opm/utility/graph/AssembledConnections.hpp>
#include <opm/utility/ECLGraph.hpp>
#include "opm/utility/ECLWellSolution.hpp"
#include "opm/flowdiagnostics/ConnectivityGraph.hpp"
#include "opm/flowdiagnostics/ConnectionValues.hpp"
#include "opm/flowdiagnostics/CellSetValues.hpp"
#include "opm/flowdiagnostics/Toolbox.hpp"
const std::string casePath = "\\\\csfiles\\Store\\ProjectData\\StatoilReservoir\\ReferenceCases\\simple_FlowDiag_Model\\";
Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G, const int step)
{
using ConnVals = Opm::FlowDiagnostics::ConnectionValues;
using NConn = ConnVals::NumConnections;
using NPhas = ConnVals::NumPhases;
const auto nconn = NConn{ G.numConnections() };
const auto nphas = NPhas{ 3 };
auto flux = ConnVals(nconn, nphas);
auto phas = ConnVals::PhaseID{ 0 };
for(const auto& p :{ Opm::BlackoilPhases::Aqua ,
Opm::BlackoilPhases::Liquid ,
Opm::BlackoilPhases::Vapour })
{
const auto pflux = G.flux(p, step);
if(! pflux.empty())
{
assert (pflux.size() == nconn.total);
auto conn = ConnVals::ConnID{ 0 };
for(const auto& v : pflux)
{
flux(conn, phas) = v;
conn.id += 1;
}
}
phas.id += 1;
}
return flux;
}
Opm::FlowDiagnostics::Toolbox
initialiseFlowDiagnostics(const Opm::ECLGraph& G,
const std::vector<Opm::ECLWellSolution::WellData>& well_fluxes,
const int step)
{
const auto connGraph = Opm::FlowDiagnostics::
ConnectivityGraph{ static_cast<int>(G.numCells()),
G.neighbours() };
using FDT = Opm::FlowDiagnostics::Toolbox;
auto fl = extractFluxField(G, step);
const size_t num_conn = fl.numConnections();
const size_t num_phases = fl.numPhases();
for(size_t conn = 0; conn < num_conn; ++conn)
{
using Co = Opm::FlowDiagnostics::ConnectionValues::ConnID;
using Ph = Opm::FlowDiagnostics::ConnectionValues::PhaseID;
for(size_t phase = 0; phase < num_phases; ++phase)
{
fl(Co{ conn }, Ph{ phase }) /= 86400; // HACK! converting to SI.
}
}
Opm::FlowDiagnostics::CellSetValues inflow;
for(const auto& well : well_fluxes)
{
for(const auto& completion : well.completions)
{
const int grid_index = completion.grid_index;
const auto& ijk = completion.ijk;
const int cell_index = G.activeCell(ijk, grid_index);
inflow.addCellValue(cell_index, completion.reservoir_inflow_rate);
}
}
// Create the Toolbox.
auto tool = FDT{ connGraph };
tool.assignPoreVolume(G.poreVolume());
tool.assignConnectionFlux(fl);
tool.assignInflowFlux(inflow);
return tool;
}
#include "exampleSetup.hpp"
TEST(opm_flowdiagnostics_test, basic_construction)
{
auto g = Opm::AssembledConnections{};
auto s = Opm::FlowDiagnostics::CellSet{};
try
try
{
Opm::ECLGraph graph = Opm::ECLGraph::load(casePath + "SIMPLE.EGRID",
casePath + "SIMPLE.INIT");
graph.assignFluxDataSource(casePath + "SIMPLE.UNRST");
int step = 2;
if ( ! graph.selectReportStep(step) )
{
std::ostringstream os;
Opm::ECLGraph eclGraph = Opm::ECLGraph::load(casePath + "SIMPLE.EGRID",
casePath + "SIMPLE.INIT");
eclGraph.assignFluxDataSource(casePath + "SIMPLE.UNRST");
os << "Report Step " << step
<< " is Not Available in Result Set '"
<< casePath << '\'';
Opm::ECLWellSolution wsol(casePath + "SIMPLE.UNRST");
auto well_fluxes = wsol.solution(2, eclGraph.numGrids());
throw std::domain_error(os.str());
}
// Opm::FlowDiagnostics::ConnectivityGraph connGraph( static_cast<int>(eclGraph.numCells()),
// eclGraph.neighbours() );
auto fdTool = initialiseFlowDiagnostics(eclGraph, well_fluxes, 2);
Opm::FlowDiagnostics::Toolbox fdTool = example::initialiseFlowDiagnostics(graph);
// Solve for time of flight.
using FDT = Opm::FlowDiagnostics::Toolbox;
std::vector<Opm::FlowDiagnostics::CellSet> start;
auto sol = fdTool.computeInjectionDiagnostics(start);
std::vector<double> globalTimeOfFlight = sol.fd.timeOfFlight();
const auto& tof = sol.fd.timeOfFlight();
// Write it to standard out.
std::cout.precision(16);
for ( double t : tof )
{
std::cout << t << '\n';
}
}
catch(const std::exception& e)
catch ( const std::exception& e )
{
std::cerr << "Caught exception: " << e.what() << '\n';
}

View File

@ -1,78 +0,0 @@
/*
Copyright 2010, 2011, 2012 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_BLACKOILPHASES_HEADER_INCLUDED
#define OPM_BLACKOILPHASES_HEADER_INCLUDED
namespace Opm
{
class BlackoilPhases
{
public:
static const int MaxNumPhases = 3;
// enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
};
struct PhaseUsage : public BlackoilPhases
{
int num_phases;
int phase_used[MaxNumPhases];
int phase_pos[MaxNumPhases];
};
/// Check or assign presence of a formed, free phase. Limited to
/// the 'BlackoilPhases'.
///
/// Use a std::vector<PhasePresence> to represent the conditions
/// in an entire model.
class PhasePresence
{
public:
PhasePresence()
: present_(0)
{}
bool hasFreeWater() const { return present(BlackoilPhases::Aqua ); }
bool hasFreeOil () const { return present(BlackoilPhases::Liquid); }
bool hasFreeGas () const { return present(BlackoilPhases::Vapour); }
void setFreeWater() { insert(BlackoilPhases::Aqua ); }
void setFreeOil () { insert(BlackoilPhases::Liquid); }
void setFreeGas () { insert(BlackoilPhases::Vapour); }
bool operator==(const PhasePresence& other) const { return present_ == other.present_; }
bool operator!=(const PhasePresence& other) const { return !this->operator==(other); }
private:
unsigned char present_;
bool present(const BlackoilPhases::PhaseIndex i) const
{
return present_ & (1 << i);
}
void insert(const BlackoilPhases::PhaseIndex i)
{
present_ |= (1 << i);
}
};
} // namespace Opm
#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED