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

Merge from upstream.

Browse Source
This commit is contained in:
Bård Skaflestad 2012-03-06 21:41:33 +01:00
commit f54a6c0979
7 changed files with 824 additions and 213 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -54,6 +54,7 @@ opm/core/utility/StopWatch.cpp \
opm/core/utility/newwells.c \
opm/core/utility/writeVtkData.cpp \
opm/core/GridManager.cpp \
opm/core/WellsManager.cpp \
opm/core/linalg/sparse_sys.c \
opm/core/linalg/LinearSolverInterface.cpp \
opm/core/pressure/cfsh.c \
@ -158,6 +159,7 @@ opm/core/well.h \
opm/core/newwells.h \
opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/WellsManager.hpp \
opm/core/pressure/fsh.h \
opm/core/pressure/HybridPressureSolver.hpp \
opm/core/pressure/IncompTpfa.hpp \

3
examples/Makefile.am
View File

@ -15,5 +15,6 @@ spu_2p_SOURCES = spu_2p.cpp
spu_2p_LDADD = \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB) \
$(LDADD)
$(LDADD) \
$(LAPACK_LIBS)
endif

439
examples/spu_2p.cpp
View File

@ -31,7 +31,7 @@
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"
@ -59,7 +59,7 @@
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/transport_source.h>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/NormSupport.hpp>
@ -96,14 +96,14 @@
class ReservoirState {
public:
ReservoirState(const UnstructuredGrid* g, const int num_phases = 2)
: press_ (g->number_of_cells, 0.0),
fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0),
sat_ (num_phases * g->number_of_cells, 0.0)
: press_ (g->number_of_cells, 0.0),
fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0),
sat_ (num_phases * g->number_of_cells, 0.0)
{
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[num_phases*cell + num_phases - 1] = 1.0;
}
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[num_phases*cell + num_phases - 1] = 1.0;
}
}
int numPhases() const { return sat_.size()/press_.size(); }
@ -129,16 +129,16 @@ private:
template <class State>
void outputState(const UnstructuredGrid* grid,
const State& state,
const int step,
const std::string& output_dir)
const State& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
@ -150,14 +150,14 @@ void outputState(const UnstructuredGrid* grid,
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
@ -168,70 +168,70 @@ class SimpleFluid2pWrappingProps
{
public:
SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props)
: props_(props),
smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases())
: props_(props),
smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases())
{
if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases.");
}
const int num_cells = props.numCells();
std::vector<int> cells(num_cells);
for (int c = 0; c < num_cells; ++c) {
cells[c] = c;
}
props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]);
if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases.");
}
const int num_cells = props.numCells();
std::vector<int> cells(num_cells);
for (int c = 0; c < num_cells; ++c) {
cells[c] = c;
}
props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]);
}
double density(int phase) const
{
return props_.density()[phase];
return props_.density()[phase];
}
template <class Sat,
class Mob,
class DMob>
class Mob,
class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity();
mob[0] /= mu[0];
mob[1] /= mu[1];
// Recall that we use Fortran ordering for kr derivatives,
// therefore dmob[i*2 + j] is row j and column i of the
// matrix.
// Each row corresponds to a kr function, so which mu to
// divide by also depends on the row, j.
dmob[0*2 + 0] /= mu[0];
dmob[0*2 + 1] /= mu[1];
dmob[1*2 + 0] /= mu[0];
dmob[1*2 + 1] /= mu[1];
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity();
mob[0] /= mu[0];
mob[1] /= mu[1];
// Recall that we use Fortran ordering for kr derivatives,
// therefore dmob[i*2 + j] is row j and column i of the
// matrix.
// Each row corresponds to a kr function, so which mu to
// divide by also depends on the row, j.
dmob[0*2 + 0] /= mu[0];
dmob[0*2 + 1] /= mu[1];
dmob[1*2 + 0] /= mu[0];
dmob[1*2 + 1] /= mu[1];
}
template <class Sat,
class Pcap,
class DPcap>
class Pcap,
class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const
{
double pc[2];
double dpc[4];
props_.capPress(1, &s[0], &c, pc, dpc);
pcap = pc[0];
ASSERT(pc[1] == 0.0);
dpcap = dpc[0];
ASSERT(dpc[1] == 0.0);
ASSERT(dpc[2] == 0.0);
ASSERT(dpc[3] == 0.0);
double pc[2];
double dpc[4];
props_.capPress(1, &s[0], &c, pc, dpc);
pcap = pc[0];
ASSERT(pc[1] == 0.0);
dpcap = dpc[0];
ASSERT(dpc[1] == 0.0);
ASSERT(dpc[2] == 0.0);
ASSERT(dpc[3] == 0.0);
}
double s_min(int c) const
{
return smin_[c];
return smin_[c];
}
double s_max(int c) const
{
return smax_[c];
return smax_[c];
}
private:
@ -258,12 +258,12 @@ public:
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
@ -284,10 +284,10 @@ main(int argc, char** argv)
const bool output = param.getDefault("output", true);
std::string output_dir;
if (output) {
output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir);
create_directories(fpath);
output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir);
create_directories(fpath);
}
// If we have a "deck_filename", grid and props will be read from that.
@ -295,25 +295,25 @@ main(int argc, char** argv)
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
// Grid init
grid.reset(new Opm::GridManager(deck));
// Rock and fluid init
const int* gc = grid->c_grid()->global_cell;
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
// Grid init
grid.reset(new Opm::GridManager(deck));
// Rock and fluid init
const int* gc = grid->c_grid()->global_cell;
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
}
// Extra rock init.
@ -329,18 +329,18 @@ main(int argc, char** argv)
double g = param.getDefault("gravity", 0.0);
bool use_gravity = g != 0.0;
if (use_gravity) {
gravity[grid->c_grid()->dimensions - 1] = g;
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
gravity[grid->c_grid()->dimensions - 1] = g;
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
bool use_segregation_split = false;
bool use_column_solver = false;
if (use_gravity && use_reorder) {
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
if (use_segregation_split) {
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
}
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
if (use_segregation_split) {
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
}
}
// Solvers init.
@ -366,7 +366,7 @@ main(int argc, char** argv)
typedef std::map<int, std::vector<int> > ColMap;
ColMap columns;
if (use_column_solver) {
Opm::extractColumn(*grid->c_grid(), columns);
Opm::extractColumn(*grid->c_grid(), columns);
}
Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid());
@ -382,57 +382,80 @@ main(int argc, char** argv)
int scenario = param.getDefault("scenario", 0);
switch (scenario) {
case 0:
{
std::cout << "==== Scenario 0: single-cell source and sink.\n";
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
break;
}
{
std::cout << "==== Scenario 0: single-cell source and sink.\n";
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
break;
}
case 1:
{
std::cout << "==== Scenario 1: half source, half sink.\n";
double flow_per_sec = 0.1*porevol[0]/Opm::unit::day;
std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec);
std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec);
break;
}
{
std::cout << "==== Scenario 1: half source, half sink.\n";
double flow_per_sec = 0.1*porevol[0]/Opm::unit::day;
std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec);
std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec);
break;
}
case 2:
{
std::cout << "==== Scenario 2: gravity convection.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
<< std::endl;
}
std::vector<double>& sat = state.saturation();
const int *glob_cell = grid->c_grid()->global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
sat[2*cell] = left ? 1.0 : 0.0;
sat[2*cell + 1] = 1.0 - sat[2*cell];
}
break;
}
{
std::cout << "==== Scenario 2: gravity convection.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
<< std::endl;
}
std::vector<double>& sat = state.saturation();
const int *glob_cell = grid->c_grid()->global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
sat[2*cell] = left ? 1.0 : 0.0;
sat[2*cell + 1] = 1.0 - sat[2*cell];
}
break;
}
case 3:
{
std::cout << "==== Scenario 2: gravity convection.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
<< std::endl;
}
std::vector<double>& sat = state.saturation();
const int *glob_cell = grid->c_grid()->global_cell;
// Heavy on top
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool top = (gc / cd[0] / cd[1]) < cd[2]/2;
sat[2*cell] = top ? 1.0 : 0.0;
sat[2*cell + 1 ] = 1.0 - sat[2*cell];
}
break;
}
default:
{
THROW("==== Scenario " << scenario << " is unknown.");
}
{
THROW("==== Scenario " << scenario << " is unknown.");
}
}
TransportSource* tsrc = create_transport_source(2, 2);
double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 };
double zdummy[] = { 0.0, 0.0 };
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc);
} else if (src[cell] < 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc);
}
if (src[cell] > 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc);
} else if (src[cell] < 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc);
}
}
std::vector<double> reorder_src = src;
@ -442,9 +465,9 @@ main(int argc, char** argv)
double current_time = 0.0;
double total_time = stepsize*num_psteps;
if (!use_reorder || use_segregation_split) {
ctrl.max_it = param.getDefault("max_it", 20);
ctrl.verbosity = param.getDefault("verbosity", 0);
ctrl.max_it_ls = param.getDefault("max_it_ls", 5);
ctrl.max_it = param.getDefault("max_it", 20);
ctrl.verbosity = param.getDefault("verbosity", 0);
ctrl.max_it_ls = param.getDefault("max_it_ls", 5);
}
// Linear solver init.
@ -455,19 +478,19 @@ main(int argc, char** argv)
// and computeTotalMobilityOmega().
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
allcells[cell] = cell;
}
// Warn if any parameters are unused.
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
// Write parameters used for later reference.
if (output) {
param.writeParam(output_dir + "/spu_2p.param");
param.writeParam(output_dir + "/spu_2p.param");
}
// Main simulation loop.
@ -480,74 +503,74 @@ main(int argc, char** argv)
std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl;
for (int pstep = 0; pstep < num_psteps; ++pstep) {
std::cout << "\n\n--------------- Simulation step number " << pstep
<< " ---------------"
<< "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
<< "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
<< "\n Total time (days) " << Opm::unit::convert::to(total_time, Opm::unit::day)
<< "\n" << std::endl;
<< " ---------------"
<< "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
<< "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
<< "\n Total time (days) " << Opm::unit::convert::to(total_time, Opm::unit::day)
<< "\n" << std::endl;
if (output) {
outputState(grid->c_grid(), state, pstep, output_dir);
}
if (output) {
outputState(grid->c_grid(), state, pstep, output_dir);
}
// Solve pressure.
if (use_gravity) {
computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega);
} else {
computeTotalMobility(*props, allcells, state.saturation(), totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
// Solve pressure.
if (use_gravity) {
computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega);
} else {
computeTotalMobility(*props, allcells, state.saturation(), totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
// Solve transport
transport_timer.start();
if (use_reorder) {
// We must treat reorder_src here,
// if we are to handle anything but simple water
// injection, since it is expected to be
// equal to total outflow (if negative)
// and water inflow (if positive).
// Also, for anything but noflow boundaries,
// boundary flows must be accumulated into
// source term following the same convention.
Opm::toWaterSat(state.saturation(), reorder_sat);
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]);
Opm::toBothSat(reorder_sat, state.saturation());
if (use_segregation_split) {
if (use_column_solver) {
colsolver.solve(columns, stepsize, state.saturation());
} else {
std::vector<double> fluxes = state.faceflux();
std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0);
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
state.faceflux() = fluxes;
}
}
} else {
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Solve transport
transport_timer.start();
if (use_reorder) {
// We must treat reorder_src here,
// if we are to handle anything but simple water
// injection, since it is expected to be
// equal to total outflow (if negative)
// and water inflow (if positive).
// Also, for anything but noflow boundaries,
// boundary flows must be accumulated into
// source term following the same convention.
Opm::toWaterSat(state.saturation(), reorder_sat);
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]);
Opm::toBothSat(reorder_sat, state.saturation());
if (use_segregation_split) {
if (use_column_solver) {
colsolver.solve(columns, stepsize, state.saturation());
} else {
std::vector<double> fluxes = state.faceflux();
std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0);
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
state.faceflux() = fluxes;
}
}
} else {
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
current_time += stepsize;
current_time += stepsize;
}
total_timer.stop();
std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
if (output) {
outputState(grid->c_grid(), state, num_psteps, output_dir);
outputState(grid->c_grid(), state, num_psteps, output_dir);
}
destroy_transport_source(tsrc);

517
opm/core/WellsManager.cpp Normal file
View File

@ -0,0 +1,517 @@
/*
Copyright 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/>.
*/
#include <opm/core/WellsManager.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <tr1/array>
#include <cmath>
// Helper structs and functions for the implementation.
namespace
{
struct WellData
{
well_type type;
control_type control;
double target;
double reference_bhp_depth;
surface_component injected_phase;
};
struct PerfData
{
int cell;
double well_index;
};
int prod_control_mode(const std::string& control)
{
const int num_prod_control_modes = 8;
static std::string prod_control_modes[num_prod_control_modes] =
{std::string("ORAT"), std::string("WRAT"), std::string("GRAT"),
std::string("LRAT"), std::string("RESV"), std::string("BHP"),
std::string("THP"), std::string("GRUP") };
int m = -1;
for (int i=0; i<num_prod_control_modes; ++i) {
if (control == prod_control_modes[i]) {
m = i;
break;
}
}
if (m >= 0) {
return m;
} else {
THROW("Unknown well control mode = " << control << " in input file");
}
}
int inje_control_mode(const std::string& control)
{
const int num_inje_control_modes = 5;
static std::string inje_control_modes[num_inje_control_modes] =
{std::string("RATE"), std::string("RESV"), std::string("BHP"),
std::string("THP"), std::string("GRUP") };
int m = -1;
for (int i=0; i<num_inje_control_modes; ++i) {
if (control == inje_control_modes[i]) {
m = i;
break;
}
}
if (m >= 0) {
return m;
} else {
THROW("Unknown well control mode = " << control << " in input file");
}
}
std::tr1::array<double, 3> getCubeDim(const UnstructuredGrid& grid, int cell)
{
using namespace std;
tr1::array<double, 3> cube;
int num_local_faces = grid.cell_facepos[cell + 1] - grid.cell_facepos[cell];
vector<double> x(num_local_faces);
vector<double> y(num_local_faces);
vector<double> z(num_local_faces);
for (int lf=0; lf<num_local_faces; ++ lf) {
int face = grid.cell_faces[grid.cell_facepos[cell] + lf];
const double* centroid = &grid.face_centroids[grid.dimensions*face];
x[lf] = centroid[0];
y[lf] = centroid[1];
z[lf] = centroid[2];
}
cube[0] = *max_element(x.begin(), x.end()) - *min_element(x.begin(), x.end());
cube[1] = *max_element(y.begin(), y.end()) - *min_element(y.begin(), y.end());
cube[2] = *max_element(z.begin(), z.end()) - *min_element(z.begin(), z.end());
return cube;
}
// Use the Peaceman well model to compute well indices.
// radius is the radius of the well.
// cubical contains [dx, dy, dz] of the cell.
// (Note that the well model asumes that each cell is a cuboid).
// cell_permeability is the permeability tensor of the given cell.
// returns the well index of the cell.
double computeWellIndex(const double radius,
const std::tr1::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor)
{
using namespace std;
// sse: Using the Peaceman model.
// NOTE: The formula is valid for cartesian grids, so the result can be a bit
// (in worst case: there is no upper bound for the error) off the mark.
const double permx = cell_permeability[0];
const double permy = cell_permeability[3*1 + 1];
double effective_perm = sqrt(permx*permy);
// sse: The formula for r_0 can be found on page 39 of
// "Well Models for Mimetic Finite Differerence Methods and Improved Representation
// of Wells in Multiscale Methods" by Ingeborg Skjelkvåle Ligaarden.
assert(permx > 0.0);
assert(permy > 0.0);
double kxoy = permx / permy;
double kyox = permy / permx;
double r0_denominator = pow(kyox, 0.25) + pow(kxoy, 0.25);
double r0_numerator = sqrt((sqrt(kyox)*cubical[0]*cubical[0]) +
(sqrt(kxoy)*cubical[1]*cubical[1]));
assert(r0_denominator > 0.0);
double r0 = 0.28 * r0_numerator / r0_denominator;
assert(radius > 0.0);
assert(r0 > 0.0);
if (r0 < radius) {
std::cout << "ERROR: Too big well radius detected.";
std::cout << "Specified well radius is " << radius
<< " while r0 is " << r0 << ".\n";
}
const long double two_pi = 6.2831853071795864769252867665590057683943387987502116419498;
double wi_denominator = log(r0 / radius) + skin_factor;
double wi_numerator = two_pi * cubical[2];
assert(wi_denominator > 0.0);
double wi = effective_perm * wi_numerator / wi_denominator;
assert(wi > 0.0);
return wi;
}
} // anonymous namespace
namespace Opm
{
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability)
: w_(0)
{
if (grid.dimensions != 3) {
THROW("We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
}
// NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells.
std::vector<std::string> keywords;
keywords.push_back("WELSPECS");
keywords.push_back("COMPDAT");
// keywords.push_back("WELTARG");
if (!deck.hasFields(keywords)) {
MESSAGE("Missing well keywords in deck, initializing no wells.");
return;
}
if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) {
THROW("Needed field is missing in file");
}
// These data structures will be filled in this constructor,
// then used to initialize the Wells struct.
std::vector<std::string> well_names;
std::vector<WellData> well_data;
std::vector<std::vector<PerfData> > wellperf_data;
// Get WELSPECS data
const WELSPECS& welspecs = deck.getWELSPECS();
const int num_wells = welspecs.welspecs.size();
well_names.reserve(num_wells);
well_data.reserve(num_wells);
wellperf_data.resize(num_wells);
for (int w = 0; w < num_wells; ++w) {
well_names.push_back(welspecs.welspecs[w].name_);
WellData wd;
well_data.push_back(wd);
well_data.back().reference_bhp_depth = welspecs.welspecs[w].datum_depth_BHP_;
if (welspecs.welspecs[w].datum_depth_BHP_ < 0.0) {
// Set refdepth to a marker value, will be changed
// after getting perforation data to the centroid of
// the cell of the top well perforation.
well_data.back().reference_bhp_depth = -1e100;
}
}
// global_cell is a map from compressed cells to Cartesian grid cells.
// We must make the inverse lookup.
const int* global_cell = grid.global_cell;
const int* cpgdim = grid.cartdims;
std::map<int,int> cartesian_to_compressed;
for (int i = 0; i < grid.number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
// Get COMPDAT data
const COMPDAT& compdat = deck.getCOMPDAT();
const int num_compdat = compdat.compdat.size();
for (int kw = 0; kw < num_compdat; ++kw) {
// Extract well name, or the part of the well name that
// comes before the '*'.
std::string name = compdat.compdat[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
// Look for well with matching name.
bool found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { // equal
// We have a matching name.
int ix = compdat.compdat[kw].grid_ind_[0] - 1;
int jy = compdat.compdat[kw].grid_ind_[1] - 1;
int kz1 = compdat.compdat[kw].grid_ind_[2] - 1;
int kz2 = compdat.compdat[kw].grid_ind_[3] - 1;
for (int kz = kz1; kz <= kz2; ++kz) {
int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz);
std::map<int, int>::const_iterator cgit =
cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
THROW("Cell with i,j,k indices " << ix << ' ' << jy << ' '
<< kz << " not found in grid!");
}
int cell = cgit->second;
PerfData pd;
pd.cell = cell;
if (compdat.compdat[kw].connect_trans_fac_ > 0.0) {
pd.well_index = compdat.compdat[kw].connect_trans_fac_;
} else {
double radius = 0.5*compdat.compdat[kw].diameter_;
if (radius <= 0.0) {
radius = 0.5*unit::feet;
MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
std::tr1::array<double, 3> cubical = getCubeDim(grid, cell);
const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell];
pd.well_index = computeWellIndex(radius, cubical, cell_perm,
compdat.compdat[kw].skin_factor_);
}
wellperf_data[wix].push_back(pd);
}
found = true;
break;
}
}
if (!found) {
THROW("Undefined well name: " << compdat.compdat[kw].well_
<< " in COMPDAT");
}
}
// Set up reference depths that were defaulted. Count perfs.
int num_perfs = 0;
for (int w = 0; w < num_wells; ++w) {
num_perfs += wellperf_data[w].size();
if (well_data[w].reference_bhp_depth == -1e100) {
// It was defaulted. Set reference depth to minimum perforation depth.
double min_depth = 1e100;
int num_wperfs = wellperf_data[w].size();
for (int perf = 0; perf < num_wperfs; ++perf) {
double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2];
min_depth = std::min(min_depth, depth);
}
well_data[w].reference_bhp_depth = min_depth;
}
}
// Get WCONINJE data
if (deck.hasField("WCONINJE")) {
const WCONINJE& wconinjes = deck.getWCONINJE();
const int num_wconinjes = wconinjes.wconinje.size();
for (int kw = 0; kw < num_wconinjes; ++kw) {
// Extract well name, or the part of the well name that
// comes before the '*'.
std::string name = wconinjes.wconinje[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data[wix].type = INJECTOR;
int m = inje_control_mode(wconinjes.wconinje[kw].control_mode_);
switch(m) {
case 0: // RATE
well_data[wix].control = RATE;
well_data[wix].target = wconinjes.wconinje[kw].surface_flow_max_rate_;
break;
case 1: // RESV
well_data[wix].control = RATE;
well_data[wix].target = wconinjes.wconinje[kw].fluid_volume_max_rate_;
break;
case 2: // BHP
well_data[wix].control = BHP;
well_data[wix].target = wconinjes.wconinje[kw].BHP_limit_;
break;
case 3: // THP
well_data[wix].control = BHP;
well_data[wix].target = wconinjes.wconinje[kw].THP_limit_;
break;
default:
THROW("Unknown well control mode; WCONIJE = "
<< wconinjes.wconinje[kw].control_mode_
<< " in input file");
}
if (wconinjes.wconinje[kw].injector_type_ == "WATER") {
well_data[wix].injected_phase = WATER;
} else if (wconinjes.wconinje[kw].injector_type_ == "OIL") {
well_data[wix].injected_phase = OIL;
} else if (wconinjes.wconinje[kw].injector_type_ == "GAS") {
well_data[wix].injected_phase = GAS;
} else {
THROW("Error in injector specification, found no known fluid type.");
}
}
}
if (!well_found) {
THROW("Undefined well name: " << wconinjes.wconinje[kw].well_
<< " in WCONINJE");
}
}
}
// Get WCONPROD data
if (deck.hasField("WCONPROD")) {
const WCONPROD& wconprods = deck.getWCONPROD();
const int num_wconprods = wconprods.wconprod.size();
for (int kw = 0; kw < num_wconprods; ++kw) {
std::string name = wconprods.wconprod[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data[wix].type = PRODUCER;
int m = prod_control_mode(wconprods.wconprod[kw].control_mode_);
switch(m) {
case 0: // ORAT
well_data[wix].control = RATE;
well_data[wix].target = wconprods.wconprod[kw].oil_max_rate_;
break;
case 1: // WRAT
well_data[wix].control = RATE;
well_data[wix].target = wconprods.wconprod[kw].water_max_rate_;
break;
case 2: // GRAT
well_data[wix].control = RATE;
well_data[wix].target = wconprods.wconprod[kw].gas_max_rate_;
break;
case 3: // LRAT
well_data[wix].control = RATE;
well_data[wix].target = wconprods.wconprod[kw].liquid_max_rate_;
break;
case 4: // RESV
well_data[wix].control = RATE;
well_data[wix].target = wconprods.wconprod[kw].fluid_volume_max_rate_;
break;
case 5: // BHP
well_data[wix].control = BHP;
well_data[wix].target = wconprods.wconprod[kw].BHP_limit_;
break;
case 6: // THP
well_data[wix].control = BHP;
well_data[wix].target = wconprods.wconprod[kw].THP_limit_;
break;
default:
THROW("Unknown well control mode; WCONPROD = "
<< wconprods.wconprod[kw].control_mode_
<< " in input file");
}
}
}
if (!well_found) {
THROW("Undefined well name: " << wconprods.wconprod[kw].well_
<< " in WCONPROD");
}
}
}
// Get WELTARG data
if (deck.hasField("WELTARG")) {
const WELTARG& weltargs = deck.getWELTARG();
const int num_weltargs = weltargs.weltarg.size();
for (int kw = 0; kw < num_weltargs; ++kw) {
std::string name = weltargs.weltarg[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data[wix].target = weltargs.weltarg[kw].new_value_;
break;
}
}
if (!well_found) {
THROW("Undefined well name: " << weltargs.weltarg[kw].well_
<< " in WELTARG");
}
}
}
// Debug output.
#define EXTRA_OUTPUT
#ifdef EXTRA_OUTPUT
std::cout << "\t WELL DATA" << std::endl;
for(int i = 0; i< num_wells; ++i) {
std::cout << i << ": " << well_data[i].type << " "
<< well_data[i].control << " " << well_data[i].target
<< std::endl;
}
std::cout << "\n\t PERF DATA" << std::endl;
for(int i=0; i< int(wellperf_data.size()); ++i) {
for(int j=0; j< int(wellperf_data[i].size()); ++j) {
std::cout << i << ": " << wellperf_data[i][j].cell << " "
<< wellperf_data[i][j].well_index << std::endl;
}
}
#endif
// Set up the Wells struct.
w_ = wells_create(num_wells, num_perfs);
if (!w_) {
THROW("Failed creating Wells struct.");
}
double fracs[3][3] = { { 1.0, 0.0, 0.0 },
{ 0.0, 1.0, 0.0 },
{ 0.0, 0.0, 1.0 } };
for (int w = 0; w < num_wells; ++w) {
int nperf = wellperf_data[w].size();
std::vector<int> cells(nperf);
std::vector<double> wi(nperf);
for (int perf = 0; perf < nperf; ++perf) {
cells[perf] = wellperf_data[w][perf].cell;
wi[perf] = wellperf_data[w][perf].well_index;
}
const double* zfrac = (well_data[w].type == INJECTOR) ? fracs[well_data[w].injected_phase] : 0;
int ok = wells_add(well_data[w].type, well_data[w].reference_bhp_depth, nperf,
zfrac, &cells[0], &wi[0], w_);
if (!ok) {
THROW("Failed to add a well.");
}
// We only append a single control at this point.
// TODO: Handle multiple controls.
ok = well_controls_append(well_data[w].control, well_data[w].target, w_->ctrls[w]);
if (!ok) {
THROW("Failed to add well controls.");
}
}
}
/// Destructor.
WellsManager::~WellsManager()
{
wells_destroy(w_);
}
/// Access the managed Wells.
/// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct.
const Wells* WellsManager::c_wells() const
{
return w_;
}
} // namespace Opm

68
opm/core/WellsManager.hpp Normal file
View File

@ -0,0 +1,68 @@
/*
Copyright 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_WELLSMANAGER_HEADER_INCLUDED
#define OPM_WELLSMANAGER_HEADER_INCLUDED
struct Wells;
struct UnstructuredGrid;
namespace Opm
{
class EclipseGridParser;
/// This class manages an Wells in the sense that it
/// encapsulates creation and destruction of the wells
/// data structure.
/// The resulting Wells is available through the c_wells() method.
class WellsManager
{
public:
/// Construct from input deck and grid.
/// The permeability argument may be zero if the input contain
/// well productivity indices, otherwise it must be given in
/// order to approximate these by the Peaceman formula.
WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability);
/// Destructor.
~WellsManager();
/// Access the managed Wells.
/// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct.
const Wells* c_wells() const;
private:
// Disable copying and assignment.
WellsManager(const WellsManager& other);
WellsManager& operator=(const WellsManager& other);
// The managed Wells.
Wells* w_;
};
} // namespace Opm
#endif // OPM_WELLSMANAGER_HEADER_INCLUDED

4
opm/core/utility/miscUtilities.cpp
View File

@ -47,7 +47,7 @@ namespace Opm
/// @brief Computes total mobility for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[int] cells cells with which the saturation values are associated
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::IncompPropertiesInterface& props,
@ -73,7 +73,7 @@ namespace Opm
/// @brief Computes total mobility and omega for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[int] cells cells with which the saturation values are associated
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobility
/// @param[out] omega mobility-weighted (or fractional-flow weighted)

4
opm/core/utility/miscUtilities.hpp
View File

@ -40,7 +40,7 @@ namespace Opm
/// @brief Computes total mobility for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[int] cells cells with which the saturation values are associated
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::IncompPropertiesInterface& props,
@ -50,7 +50,7 @@ namespace Opm
/// @brief Computes total mobility and omega for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[int] cells cells with which the saturation values are associated
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobility
/// @param[out] omega mobility-weighted (or fractional-flow weighted)