This commit is contained in:
Xavier Raynaud committed 2012-04-10 17:18:38 +02:00
commit 820b5d25a7
21 files changed
+115 -113

No files matched your search

+17 -2
View File
@@ -3,8 +3,11 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/WellsManager.hpp>
#include <opm/core/GridManager.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include "opm/core/newwells.h"
#include "opm/core/grid.h"
int main(int argc, char** argv) {
@@ -12,7 +15,7 @@ int main(int argc, char** argv) {
using namespace Opm;
ParameterGroup parameters( argc, argv, false );
std::string file_name = parameters.getDefault<std::string>("inputdeck", "data.data");
// Read input file
EclipseGridParser parser(file_name);
std::cout << "Done!" << std::endl;
@@ -22,6 +25,18 @@ int main(int argc, char** argv) {
// Finally handle the wells
WellsManager wells(parser, *grid.c_grid(), NULL);
std::vector<int> global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells);
std::cout << "ahoi" << std::endl;
double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
IncompPropertiesFromDeck incomp_properties(parser, global_cells);
std::cout << "there" << std::endl;
LinearSolverUmfpack umfpack_solver;
std::cout << "here" << std::endl;
IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(),
gravity, umfpack_solver, wells.c_wells());
return 0;
}
+10 -4
View File
@@ -41,11 +41,11 @@ namespace Opm
}
std::tr1::shared_ptr<WellsGroupInterface> child;
for (int i = 0; i < roots_.size(); ++i) {
for (size_t i = 0; i < roots_.size(); ++i) {
if (roots_[i]->name() == child_name) {
child = roots_[i];
// We've found a new parent to the previously thought root, need to remove it
for(int j = i; j < roots_.size() - 1; ++j) {
for(size_t j = i; j < roots_.size() - 1; ++j) {
roots_[j] = roots_[j+1];
}
@@ -63,13 +63,19 @@ namespace Opm
}
parent_as_group->addChild(child);
if(child->isLeafNode()) {
leaf_nodes_.push_back(child);
}
}
const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& WellCollection::getLeafNodes() const {
return leaf_nodes_;
}
WellsGroupInterface* WellCollection::findNode(std::string name)
{
for (int i = 0; i < roots_.size(); i++) {
for (size_t i = 0; i < roots_.size(); i++) {
WellsGroupInterface* result = roots_[i]->findGroup(name);
if (result) {
return result;
+5
View File
@@ -38,9 +38,14 @@ namespace Opm
void addChild(std::string child, std::string parent,
const EclipseGridParser& deck);
const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& getLeafNodes() const;
private:
// To account for the possibility of a forest
std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;
// This will be used to traverse the bottom nodes.
std::vector<std::tr1::shared_ptr<WellsGroupInterface> > leaf_nodes_;
WellsGroupInterface* findNode(std::string name);
};
+14 -6
View File
@@ -31,13 +31,17 @@ namespace Opm
: WellsGroupInterface(name, prod_spec, inj_spec)
{
}
bool WellsGroupInterface::isLeafNode() const {
return false;
}
WellsGroupInterface* WellsGroup::findGroup(std::string name_of_node)
{
if (name() == name_of_node) {
return this;
} else {
for (int i = 0; i < children_.size(); i++) {
for (size_t i = 0; i < children_.size(); i++) {
WellsGroupInterface* result = children_[i]->findGroup(name_of_node);
if (result) {
return result;
@@ -68,6 +72,10 @@ namespace Opm
}
}
bool WellNode::isLeafNode() const {
return true;
}
surface_component toSurfaceComponent(std::string type)
{
@@ -171,7 +179,7 @@ namespace Opm
bool isWell = false;
if (deck.hasField("WELSPECS")) {
WELSPECS wspecs = deck.getWELSPECS();
for (int i = 0; i < wspecs.welspecs.size(); i++) {
for (size_t i = 0; i < wspecs.welspecs.size(); i++) {
if (wspecs.welspecs[i].name_ == name) {
isWell = true;
break;
@@ -185,7 +193,7 @@ namespace Opm
InjectionSpecification injection_specification;
if (deck.hasField("WCONINJE")) {
WCONINJE wconinje = deck.getWCONINJE();
for (int i = 0; i < wconinje.wconinje.size(); i++) {
for (size_t i = 0; i < wconinje.wconinje.size(); i++) {
if (wconinje.wconinje[i].well_ == name) {
WconinjeLine line = wconinje.wconinje[i];
injection_specification.BHP_limit_ = line.BHP_limit_;
@@ -199,7 +207,7 @@ namespace Opm
ProductionSpecification production_specification;
if (deck.hasField("WCONPROD")) {
WCONPROD wconprod = deck.getWCONPROD();
for (int i = 0; i < wconprod.wconprod.size(); i++) {
for (size_t i = 0; i < wconprod.wconprod.size(); i++) {
if (wconprod.wconprod[i].well_ == name) {
WconprodLine line = wconprod.wconprod[i];
production_specification.BHP_limit_ = line.BHP_limit_;
@@ -216,7 +224,7 @@ namespace Opm
InjectionSpecification injection_specification;
if (deck.hasField("GCONINJE")) {
GCONINJE gconinje = deck.getGCONINJE();
for (int i = 0; i < gconinje.gconinje.size(); i++) {
for (size_t i = 0; i < gconinje.gconinje.size(); i++) {
if (gconinje.gconinje[i].group_ == name) {
GconinjeLine line = gconinje.gconinje[i];
injection_specification.injector_type_ = toSurfaceComponent(line.injector_type_);
@@ -229,7 +237,7 @@ namespace Opm
ProductionSpecification production_specification;
if (deck.hasField("GCONPROD")) {
GCONPROD gconprod = deck.getGCONPROD();
for (int i = 0; i < gconprod.gconprod.size(); i++) {
for (size_t i = 0; i < gconprod.gconprod.size(); i++) {
if (gconprod.gconprod[i].group_ == name) {
GconprodLine line = gconprod.gconprod[i];
production_specification.oil_max_rate_ = line.oil_max_rate_;
+10
View File
@@ -17,9 +17,17 @@ namespace Opm
InjectionSpecification inj_spec);
virtual ~WellsGroupInterface();
/// The unique identifier for the well or well group.
const std::string& name();
/// Production specifications for the well or well group.
const ProductionSpecification& prodSpec() const;
/// Injection specifications for the well or well group.
const InjectionSpecification& injSpec() const;
/// \returns true if the object is a leaf node (WellNode), false otherwise.
virtual bool isLeafNode() const;
/// \returns the pointer to the WellsGroupInterface with the given name. NULL if
/// the name is not found.a
@@ -51,6 +59,8 @@ namespace Opm
InjectionSpecification inj_spec);
virtual WellsGroupInterface* findGroup(std::string name_of_node);
virtual bool isLeafNode() const;
};
+4 -4
View File
@@ -505,22 +505,22 @@ namespace Opm
}
}
WellCollection well_collection;
WellCollection wells_;
if (deck.hasField("GRUPTREE")) {
std::cout << "Found gruptree" << std::endl;
const GRUPTREE& gruptree = deck.getGRUPTREE();
std::map<std::string, std::string>::const_iterator it = gruptree.tree.begin();
for( ; it != gruptree.tree.end(); ++it) {
well_collection.addChild(it->first, it->second, deck);
wells_.addChild(it->first, it->second, deck);
}
}
if(deck.hasField("WELSPECS")) {
WELSPECS welspecs = deck.getWELSPECS();
for(int i = 0; i < welspecs.welspecs.size(); ++i) {
for(size_t i = 0; i < welspecs.welspecs.size(); ++i) {
WelspecsLine line = welspecs.welspecs[i];
well_collection.addChild(line.name_, line.group_, deck);
wells_.addChild(line.name_, line.group_, deck);
}
}
+3 -1
View File
@@ -19,7 +19,7 @@
#ifndef OPM_WELLSMANAGER_HEADER_INCLUDED
#define OPM_WELLSMANAGER_HEADER_INCLUDED
#include <opm/core/WellCollection.hpp>
struct Wells;
struct UnstructuredGrid;
@@ -61,6 +61,8 @@ namespace Opm
WellsManager& operator=(const WellsManager& other);
// The managed Wells.
Wells* w_;
WellCollection well_collection_;
};
+3 -57
View File
@@ -125,7 +125,7 @@ namespace Opm
double B[2]; // Must be enough since component classes do not handle more than 2.
pvt_.B(1, 0, 0, B);
// Compute A matrix
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = A + i*np*np;
std::fill(m, m + np*np, 0.0);
@@ -137,66 +137,12 @@ namespace Opm
// Derivative of A matrix.
if (dAdp) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = dAdp + i*np*np;
std::fill(m, m + np*np, 0.0);
}
}
#if 0
// This is copied from BlackoilPropertiesFromDeck.
const int np = numPhases();
B_.resize(n*np);
R_.resize(n*np);
if (dAdp) {
dB_.resize(n*np);
dR_.resize(n*np);
pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
} else {
pvt_.B(n, p, z, &B_[0]);
pvt_.R(n, p, z, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
pvt_.phaseUsed()[BlackoilPhases::Vapour];
const int o = phase_pos[BlackoilPhases::Liquid];
const int g = phase_pos[BlackoilPhases::Vapour];
// Compute A matrix
#pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = A + i*np*np;
std::fill(m, m + np*np, 0.0);
// Diagonal entries.
for (int phase = 0; phase < np; ++phase) {
m[phase + phase*np] = 1.0/B_[i*np + phase];
}
// Off-diagonal entries.
if (oil_and_gas) {
m[o + g*np] = R_[i*np + g]/B_[i*np + g];
m[g + o*np] = R_[i*np + o]/B_[i*np + o];
}
}
// Derivative of A matrix.
if (dAdp) {
#pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = dAdp + i*np*np;
std::fill(m, m + np*np, 0.0);
// Diagonal entries.
for (int phase = 0; phase < np; ++phase) {
m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase];
}
// Off-diagonal entries.
if (oil_and_gas) {
m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g];
m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o];
}
}
}
#endif
}
/// \param[in] n Number of data points.
@@ -211,7 +157,7 @@ namespace Opm
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
@@ -131,7 +131,7 @@ namespace Opm
const int g = phase_pos[BlackoilPhases::Vapour];
// Compute A matrix
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = A + i*np*np;
std::fill(m, m + np*np, 0.0);
@@ -148,7 +148,7 @@ namespace Opm
// Derivative of A matrix.
if (dAdp) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = dAdp + i*np*np;
std::fill(m, m + np*np, 0.0);
@@ -177,7 +177,7 @@ namespace Opm
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
+3 -3
View File
@@ -79,7 +79,7 @@ namespace Opm
{
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[np*i + phase] = viscosity_[phase];
}
@@ -93,7 +93,7 @@ namespace Opm
{
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase];
}
@@ -108,7 +108,7 @@ namespace Opm
{
const int np = numPhases();
for (int phase = 0; phase < np; ++phase) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[np*i + phase] = formation_volume_factor_[phase];
output_dBdp[np*i + phase] = 0.0;
+1 -1
View File
@@ -42,7 +42,7 @@ namespace Opm
permeability_.clear();
const int dsq = dimensions*dimensions;
permeability_.resize(num_cells*dsq, 0.0);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < num_cells; ++i) {
for (int d = 0; d < dimensions; ++d) {
permeability_[dsq*i + dimensions*d + d] = perm;
+2 -2
View File
@@ -71,13 +71,13 @@ namespace Opm
const double* s, double* kr, double* dkrds, Fun fun)
{
if (dkrds == 0) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n*np; ++i) {
kr[i] = fun.kr(s[i]);
}
return;
}
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0);
for (int phase = 0; phase < np; ++phase) {
+4 -4
View File
@@ -111,12 +111,12 @@ namespace Opm
{
const int np = phase_usage_.num_phases;
if (dkrds) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
} else {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalKr(s + np*i, kr + np*i);
}
@@ -142,12 +142,12 @@ namespace Opm
{
const int np = phase_usage_.num_phases;
if (dpcds) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
} else {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalPc(s + np*i, pc + np*i);
}
@@ -134,7 +134,7 @@ namespace Opm
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->mu(n, p, z, &data1_[0]);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
}
@@ -149,7 +149,7 @@ namespace Opm
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->B(n, p, z, &data1_[0]);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
}
@@ -166,7 +166,7 @@ namespace Opm
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
output_dBdp[phase_usage_.num_phases*i + phase] = data2_[i];
@@ -183,7 +183,7 @@ namespace Opm
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->R(n, p, z, &data1_[0]);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
}
@@ -200,7 +200,7 @@ namespace Opm
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
output_dRdp[phase_usage_.num_phases*i + phase] = data2_[i];
@@ -81,7 +81,7 @@ namespace Opm
double* output_B) const
{
if (comp_) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(p[i] - ref_press_);
@@ -100,7 +100,7 @@ namespace Opm
{
B(n, p, 0, output_B);
if (comp_) {
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_dBdp[i] = -comp_*output_B[i];
}
+3 -3
View File
@@ -75,7 +75,7 @@ namespace Opm
const double* /*z*/,
double* output_mu) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
}
@@ -86,7 +86,7 @@ namespace Opm
const double* /*z*/,
double* output_B) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = 1.0/one_over_B_(p[i]);
}
@@ -99,7 +99,7 @@ namespace Opm
double* output_dBdp) const
{
B(n, p, 0, output_B);
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(p[i]);
+5 -5
View File
@@ -91,7 +91,7 @@ namespace Opm
const double* z,
double* output_mu) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_gas(p[i], z + num_phases_*i, 2, false);
}
@@ -104,7 +104,7 @@ namespace Opm
const double* z,
double* output_B) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i);
}
@@ -119,7 +119,7 @@ namespace Opm
double* output_B,
double* output_dBdp) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
@@ -132,7 +132,7 @@ namespace Opm
const double* z,
double* output_R) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i);
}
@@ -147,7 +147,7 @@ namespace Opm
double* output_R,
double* output_dRdp) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
+5 -5
View File
@@ -168,7 +168,7 @@ namespace Opm
const double* z,
double* output_mu) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_oil(p[i], z + num_phases_*i, 2, false);
}
@@ -181,7 +181,7 @@ namespace Opm
const double* z,
double* output_B) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i);
}
@@ -196,7 +196,7 @@ namespace Opm
double* output_B,
double* output_dBdp) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
@@ -209,7 +209,7 @@ namespace Opm
const double* z,
double* output_R) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i);
}
@@ -224,7 +224,7 @@ namespace Opm
double* output_R,
double* output_dRdp) const
{
#pragma omp parallel for
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
+8 -4
View File
@@ -38,14 +38,18 @@ namespace Opm
/// and N == g.number_of_cells.
/// \param[in] gravity Gravity vector. If nonzero, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL
IncompTpfa::IncompTpfa(const UnstructuredGrid& g,
const double* permeability,
const double* gravity,
const LinearSolverInterface& linsolver)
const LinearSolverInterface& linsolver,
const struct Wells* wells)
: grid_(g),
linsolver_(linsolver),
htrans_(g.cell_facepos[ g.number_of_cells ]),
trans_ (g.number_of_faces)
trans_ (g.number_of_faces),
wells_(wells)
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_htrans_compute(gg, permeability, &htrans_[0]);
@@ -110,7 +114,7 @@ namespace Opm
}
}
ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL };
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
@@ -172,7 +176,7 @@ namespace Opm
}
}
ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL };
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
+7 -2
View File
@@ -25,6 +25,7 @@
struct UnstructuredGrid;
struct ifs_tpfa_data;
struct Wells;
struct FlowBoundaryConditions;
namespace Opm
@@ -47,10 +48,13 @@ namespace Opm
/// and N == g.number_of_cells.
/// \param[in] gravity Gravity vector. If nonzero, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL
IncompTpfa(const UnstructuredGrid& g,
const double* permeability,
const double* gravity,
const LinearSolverInterface& linsolver);
const LinearSolverInterface& linsolver,
const struct Wells* wells = 0);
/// Destructor.
~IncompTpfa();
@@ -113,7 +117,8 @@ namespace Opm
::std::vector<double> trans_ ;
::std::vector<double> gpress_;
::std::vector<double> gpress_omegaweighted_;
const struct Wells* wells_;
struct ifs_tpfa_data* h_;
};
+1
View File
@@ -22,6 +22,7 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <boost/lexical_cast.hpp>
#include <set>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <iterator>