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

Added default oil handling

Browse Source
This commit is contained in:
Kjetil Olsen Lye
2012-04-23 13:49:53 +02:00
parent fbf898b756 1f9fea92fa
commit 6670f5de17
30 changed files with 1341 additions and 346 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
Doxyfile
View File

@@ -51,7 +51,7 @@ PROJECT_LOGO =
# If a relative path is entered, it will be relative to the location
# where doxygen was started. If left blank the current directory will be used.
OUTPUT_DIRECTORY =
OUTPUT_DIRECTORY = ../Documentation
# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create
# 4096 sub-directories (in 2 levels) under the output directory of each output
@@ -610,9 +610,8 @@ WARN_LOGFILE =
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
# INPUT = tutorials/tutorial1.cpp opm/core/GridManager.hpp opm/core/utility/writeVtkData.hpp opm/core/utility/writeVtkData.cpp
INPUT = tutorials opm/core/linalg/LinearSolverUmfpack.hpp opm/core/pressure/IncompTpfa.hpp opm/core/pressure/FlowBCManager.hpp opm/core/utility/miscUtilities.hpp opm/core/GridManager.hpp opm/core/utility/writeVtkData.hpp opm/core/grid.h
#INPUT = opm/core tutorials examples
INPUT = tutorials
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is
@@ -690,7 +689,7 @@ EXAMPLE_RECURSIVE = NO
# directories that contain image that are included in the documentation (see
# the \image command).
IMAGE_PATH = Figure/
IMAGE_PATH = ../Documentation/Figure/
# The INPUT_FILTER tag can be used to specify a program that doxygen should
# invoke to filter for each input file. Doxygen will invoke the filter program

7
examples/spu_2p.cpp
View File

@@ -266,7 +266,12 @@ main(int argc, char** argv)
output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir);
create_directories(fpath);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
output_interval = param.getDefault("output_interval", output_interval);
}
const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);

78
generate_doc_figures.py
View File

@@ -1,13 +1,25 @@
from subprocess import call
from paraview.simple import *
from paraview import servermanager
from os import remove
from os import remove, mkdir, curdir
from os.path import join, isdir
figure_path = "../Documentation/Figure"
tutorial_data_path = curdir
tutorial_path = "tutorials"
collected_garbage_file = []
if not isdir(figure_path):
mkdir(figure_path)
connection = servermanager.Connect()
# tutorial 1
call("tutorials/tutorial1")
grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial1.vtu")
call(join(tutorial_path, "tutorial1"))
data_file_name = join(tutorial_data_path, "tutorial1.vtu")
grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name)
collected_garbage_file.append(data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
@@ -17,28 +29,27 @@ dp.AmbientColor = [1, 0, 0]
view = GetActiveView()
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(20, 20, 110)
camera.SetViewUp(0.5,0.3,0.7)
camera.SetPosition(4, -6, 5)
camera.SetViewUp(-0.19, 0.4, 0.9)
camera.SetViewAngle(30)
camera.SetFocalPoint(1,1,0.5)
camera.SetFocalPoint(1.5, 1.5, 1)
Render()
WriteImage("Figure/tutorial1.png")
WriteImage(join(figure_path, "tutorial1.png"))
Hide(grid)
# remove("tutorial1.vtu")
# tutorial 2
call("tutorials/tutorial2")
grid = servermanager.sources.XMLUnstructuredGridReader(FileName="tutorial2.vtu")
call(join(tutorial_path, "tutorial2"))
data_file_name = join(tutorial_data_path, "tutorial2.vtu")
grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name)
collected_garbage_file.append(data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
data = grid.GetCellDataInformation()
pressure = data.GetArray('pressure')
pmin = pressure.GetRange()[0]
pmax = pressure.GetRange()[1]
dp.LookupTable = MakeBlueToRedLT(0.5*(pmin+pmax)-0.2*(pmax-pmin), 0.5*(pmin+pmax)-0.2*(pmax+pmin))
dp.ColorArrayName = 'pressure'
pres = grid.CellData.GetArray(0)
pres_lookuptable = GetLookupTableForArray( "pressure", 1, RGBPoints=[pres.GetRange()[0], 1, 0, 0, pres.GetRange()[1], 0, 0, 1] )
dp.LookupTable = pres_lookuptable
view = GetActiveView()
view.Background = [1, 1, 1]
camera = GetActiveCamera()
@@ -47,5 +58,38 @@ camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
camera.SetFocalPoint(20, 20, 0.5)
Render()
WriteImage("Figure/tutorial2.png")
# remove("tutorial2.vtu")
WriteImage(join(figure_path, "tutorial2.png"))
Hide(grid)
# tutorial 3
call(join(tutorial_path, "tutorial3"))
for case in range(0,20):
data_file_name = join(tutorial_data_path, "tutorial3-"+"%(case)03d"%{"case": case}+".vtu")
collected_garbage_file.append(data_file_name)
cases = ["000", "005", "010", "015", "019"]
for case in cases:
data_file_name = join(tutorial_data_path, "tutorial3-"+case+".vtu")
grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
dp.ColorArrayName = 'saturation'
sat = grid.CellData.GetArray(1)
sat_lookuptable = GetLookupTableForArray( "saturation", 1, RGBPoints=[0, 1, 0, 0, 1, 0, 0, 1])
dp.LookupTable = sat_lookuptable
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(100, 100, 550)
camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
camera.SetFocalPoint(100, 100, 5)
Render()
WriteImage(join(figure_path, "tutorial3-"+case+".png"))
Hide(grid)
# remove temporary files
for f in collected_garbage_file:
remove(f)

2
opm/core/InjectionSpecification.hpp
View File

@@ -15,7 +15,7 @@ namespace Opm
InjectionSpecification();
surface_component injector_type_;
SurfaceComponent injector_type_;
ControlMode control_mode_;
double surface_flow_max_rate_;
double reinjection_fraction_target_;

4
opm/core/WellsGroup.cpp
View File

@@ -298,7 +298,7 @@ namespace Opm
namespace
{
surface_component toSurfaceComponent(std::string type)
SurfaceComponent toSurfaceComponent(std::string type)
{
if (type == "OIL") {
return OIL;
@@ -309,7 +309,7 @@ namespace Opm
if (type == "GAS") {
return GAS;
}
THROW("Unknown type " << type << ", could not convert to surface_component");
THROW("Unknown type " << type << ", could not convert to SurfaceComponent");
}
InjectionSpecification::ControlMode toInjectionControlMode(std::string type)

46
opm/core/WellsManager.cpp
View File

@@ -36,11 +36,11 @@ namespace
struct WellData
{
well_type type;
control_type control;
WellType type;
WellControlType control;
double target;
double reference_bhp_depth;
surface_component injected_phase;
SurfaceComponent injected_phase;
};
@@ -517,9 +517,19 @@ namespace Opm
// Apply guide rates:
for (size_t i = 0; i < well_data.size(); i++) {
if (well_collection_.getLeafNodes()[i]->prodSpec().control_mode_ == ProductionSpecification::GRUP) {
if (well_data[i].type == PRODUCER && (well_collection_.getLeafNodes()[i]->prodSpec().control_mode_ == ProductionSpecification::GRUP)) {
switch (well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_ ) {
case ProductionSpecification::OIL:
{
const ProductionSpecification& parent_prod_spec =
well_collection_.getLeafNodes()[i]->getParent()->prodSpec();
double guide_rate = well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_;
well_data[i].target = guide_rate*parent_prod_spec.oil_max_rate_;
well_data[i].control = RATE;
well_data[i].type = PRODUCER;
std::cout << "WARNING: Converting oil control to rate control!" << std::endl;
break;
}
case ProductionSpecification::NONE_GRT:
{
// Will use the group control type:
@@ -532,17 +542,18 @@ namespace Opm
well_data[i].control = RATE;
break;
default:
THROW("Unhandled group control mode " << parent_prod_spec.control_mode_);
THROW("Unhandled production specification control mode " << parent_prod_spec.control_mode_);
break;
}
}
default:
THROW("Unhandled production specification guide rate type "
<< well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_);
break;
}
}
if (well_collection_.getLeafNodes()[i]->injSpec().control_mode_ == InjectionSpecification::GRUP) {
if (well_data[i].type == INJECTOR && (well_collection_.getLeafNodes()[i]->injSpec().control_mode_ == InjectionSpecification::GRUP)) {
if (well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_ == ProductionSpecification::RAT) {
well_data[i].injected_phase = WATER; // Default for now.
well_data[i].control = RATE;
@@ -558,7 +569,7 @@ namespace Opm
std::cout << "Making well structs" << std::endl;
// Set up the Wells struct.
w_ = wells_create(num_wells, num_perfs);
w_ = create_wells(num_wells, num_perfs);
if (!w_) {
THROW("Failed creating Wells struct.");
}
@@ -576,32 +587,31 @@ namespace Opm
const double* zfrac = (well_data[w].type == INJECTOR) ? fracs[well_data[w].injected_phase] : 0;
// DIRTY DIRTY HACK
if(well_data[w].type == INJECTOR && well_data[w].injected_phase < 0 || well_data[w].injected_phase > 2){
if(well_data[w].type == INJECTOR && (well_data[w].injected_phase < 0 || well_data[w].injected_phase > 2)){
zfrac = fracs[WATER];
}
int ok = wells_add(well_data[w].type, well_data[w].reference_bhp_depth, nperf,
int ok = add_well(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]);
ok = append_well_controls(well_data[w].control, well_data[w].target, w_->ctrls[w]);
w_->ctrls[w]->current = 0;
if (!ok) {
THROW("Failed to add well controls.");
}
}
std::cout << "Made well struct" << std::endl;
for(size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) {
// \TODO comment this.
for (size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) {
WellNode* node = static_cast<WellNode*>(well_collection_.getLeafNodes()[i].get());
// We know that getLeafNodes() is ordered the same way as they're indexed in w_
node->setWellsPointer(w_, i);
}
}
@@ -609,7 +619,7 @@ namespace Opm
/// Destructor.
WellsManager::~WellsManager()
{
wells_destroy(w_);
destroy_wells(w_);
}

10
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@@ -27,7 +27,7 @@ namespace Opm
{
rock_.init(deck, global_cell);
pvt_.init(deck);
satprops_.init(deck);
satprops_.init(deck, global_cell);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
@@ -199,11 +199,11 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesFromDeck::relperm(const int n,
const double* s,
const int* /*cells*/,
const int* cells,
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, kr, dkrds);
satprops_.relperm(n, s, cells, kr, dkrds);
}
@@ -218,11 +218,11 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void BlackoilPropertiesFromDeck::capPress(const int n,
const double* s,
const int* /*cells*/,
const int* cells,
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, pc, dpcds);
satprops_.capPress(n, s, cells, pc, dpcds);
}

20
opm/core/fluid/IncompPropertiesBasic.cpp
View File

@@ -46,6 +46,26 @@ namespace Opm
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
IncompPropertiesBasic::IncompPropertiesBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const double por, //porosity
const double perm,
const int dim,
const int num_cells)
{
rock_.init(dim, num_cells, por, perm);
pvt_.init(num_phases, rho, mu);
satprops_.init(num_phases, relpermfunc);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesBasic::IncompPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
IncompPropertiesBasic::~IncompPropertiesBasic()
{
}

11
opm/core/fluid/IncompPropertiesBasic.hpp
View File

@@ -53,6 +53,17 @@ namespace Opm
const int dim,
const int num_cells);
/// Construct from arguments a basic two phase fluid.
IncompPropertiesBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const double porosity,
const double permeability,
const int dim,
const int num_cells);
/// Destructor.
virtual ~IncompPropertiesBasic();

16
opm/core/fluid/IncompPropertiesFromDeck.cpp
View File

@@ -31,7 +31,7 @@ namespace Opm
{
rock_.init(deck, global_cell);
pvt_.init(deck);
satprops_.init(deck);
satprops_.init(deck, global_cell);
if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
@@ -87,7 +87,7 @@ namespace Opm
/// \return Array of P density values.
const double* IncompPropertiesFromDeck::density() const
{
return pvt_.surfaceDensities();
return pvt_.reservoirDensities();
}
/// \param[in] n Number of data points.
@@ -101,11 +101,11 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::relperm(const int n,
const double* s,
const int* /*cells*/,
const int* cells,
double* kr,
double* dkrds) const
{
satprops_.relperm(n, s, kr, dkrds);
satprops_.relperm(n, s, cells, kr, dkrds);
}
@@ -120,11 +120,11 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m_01 ...)
void IncompPropertiesFromDeck::capPress(const int n,
const double* s,
const int* /*cells*/,
const int* cells,
double* pc,
double* dpcds) const
{
satprops_.capPress(n, s, pc, dpcds);
satprops_.capPress(n, s, cells, pc, dpcds);
}
@@ -136,11 +136,11 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void IncompPropertiesFromDeck::satRange(const int n,
const int* /*cells*/,
const int* cells,
double* smin,
double* smax) const
{
satprops_.satRange(n, smin, smax);
satprops_.satRange(n, cells, smin, smax);
}
} // namespace Opm

14
opm/core/fluid/PvtPropertiesBasic.cpp
View File

@@ -59,6 +59,20 @@ namespace Opm
}
}
void PvtPropertiesBasic::init(const int num_phases,
const std::vector<double>& rho,
const std::vector<double>& visc)
{
if (num_phases > 3 || num_phases < 1) {
THROW("PvtPropertiesBasic::init() illegal num_phases: " << num_phases);
}
// We currently do not allow the user to set B.
formation_volume_factor_.clear();
formation_volume_factor_.resize(num_phases, 1.0);
density_ = rho;
viscosity_ = visc;
}
const double* PvtPropertiesBasic::surfaceDensities() const
{
return &density_[0];

6
opm/core/fluid/PvtPropertiesBasic.hpp
View File

@@ -44,6 +44,12 @@ namespace Opm
/// mu1 [mu2, mu3] (1.0) Viscosity in cP
void init(const parameter::ParameterGroup& param);
/// Initialize from arguments.
/// Basic multi phase fluid pvt properties.
void init(const int num_phases,
const std::vector<double>& rho,
const std::vector<double>& visc);
/// Number of active phases.
int numPhases() const;

18
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
View File

@@ -51,18 +51,23 @@ namespace Opm
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else {
THROW("Input is missing DENSITY\n");
}
// Make reservoir densities the same as surface densities initially.
// We will modify them with formation volume factors if found.
reservoir_density_ = surface_density_;
// Water viscosity.
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
MESSAGE("Compressibility effects in PVTW are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
// Eclipse 100 default.
@@ -76,6 +81,7 @@ namespace Opm
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
MESSAGE("Compressibility effects in PVCDO are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {
THROW("Input is missing PVCDO\n");
@@ -84,7 +90,13 @@ namespace Opm
const double* PvtPropertiesIncompFromDeck::surfaceDensities() const
{
return density_.data();
return surface_density_.data();
}
const double* PvtPropertiesIncompFromDeck::reservoirDensities() const
{
return reservoir_density_.data();
}

18
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
View File

@@ -31,9 +31,6 @@ namespace Opm
/// eclipse input (keywords DENSITY, PVTW, PVCDO).
///
/// All phases are incompressible and have constant viscosities.
/// For all the methods, the following apply: p and z are unused.
/// Output arrays shall be of size n*numPhases(), and must be valid
/// before calling the method.
/// NOTE: This class is intentionally similar to BlackoilPvtProperties.
class PvtPropertiesIncompFromDeck
{
@@ -51,11 +48,24 @@ namespace Opm
/// \return Array of size numPhases().
const double* surfaceDensities() const;
/// Densities of stock components at reservoir conditions.
/// Note: a reasonable question to ask is why there can be
/// different densities at surface and reservoir conditions,
/// when the phases are assumed incompressible. The answer is
/// that even if we approximate the phases as being
/// incompressible during simulation, the density difference
/// between surface and reservoir may be larger. For accurate
/// reporting and using data given in terms of surface values,
/// we need to handle this difference.
/// \return Array of size numPhases().
const double* reservoirDensities() const;
/// Viscosities.
const double* viscosity() const;
private:
std::tr1::array<double, 2> density_;
std::tr1::array<double, 2> surface_density_;
std::tr1::array<double, 2> reservoir_density_;
std::tr1::array<double, 2> viscosity_;
};

13
opm/core/fluid/SaturationPropsBasic.hpp
View File

@@ -45,6 +45,16 @@ namespace Opm
/// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic".
void init(const parameter::ParameterGroup& param);
enum RelPermFunc { Constant, Linear, Quadratic };
/// Initialize from arguments a basic Saturation property.
void init(const int num_phases,
const RelPermFunc& relperm_func)
{
num_phases_ = num_phases;
relperm_func_ = relperm_func;
}
/// \return P, the number of phases.
int numPhases() const;
@@ -83,8 +93,9 @@ namespace Opm
void satRange(const int n,
double* smin,
double* smax) const;
private:
enum RelPermFunc { Constant, Linear, Quadratic };
int num_phases_;
RelPermFunc relperm_func_;
};

210
opm/core/fluid/SaturationPropsFromDeck.cpp
View File

@@ -32,7 +32,8 @@ namespace Opm
}
/// Initialize from deck.
void SaturationPropsFromDeck::init(const EclipseGridParser& deck)
void SaturationPropsFromDeck::init(const EclipseGridParser& deck,
const std::vector<int>& global_cell)
{
phase_usage_ = phaseUsageFromDeck(deck);
@@ -41,49 +42,48 @@ namespace Opm
if (!phase_usage_.phase_used[Liquid]) {
THROW("SaturationPropsFromDeck::init() -- oil phase must be active.");
}
const int samples = 200;
double swco = 0.0;
double swmax = 1.0;
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
int num_cells = global_cell.size();
cell_to_func_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cell_to_func_[cell] = satnum[global_cell[cell]] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
if (swof_table.size() != 1) {
THROW("We must have exactly one SWOF table.");
num_tables = swof_table.size();
if (num_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
const std::vector<double>& sw = swof_table[0][0];
const std::vector<double>& krw = swof_table[0][1];
const std::vector<double>& krow = swof_table[0][2];
const std::vector<double>& pcow = swof_table[0][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage_.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage_.phase_pos[Aqua]] = sw.back();
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
if (sgof_table.size() != 1) {
THROW("We must have exactly one SGOF table.");
int num_sgof_tables = sgof_table.size();
if (num_sgof_tables < satfuncs_expected) {
THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
THROW("Inconsistent number of tables in SWOF and SGOF.");
}
const std::vector<double>& sg = sgof_table[0][0];
const std::vector<double>& krg = sgof_table[0][1];
const std::vector<double>& krog = sgof_table[0][2];
const std::vector<double>& pcog = sgof_table[0][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
smin_[phase_usage_.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-2) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage_.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage_.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage_.phase_pos[Liquid]] = 1.0 - swco;
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_);
}
}
@@ -101,6 +101,7 @@ namespace Opm
/// Relative permeability.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of nP relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of nP^2 relperm derivative values,
/// array must be valid before calling.
@@ -109,6 +110,7 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsFromDeck::relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const
{
@@ -116,12 +118,12 @@ namespace Opm
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalKr(s + np*i, kr + np*i);
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
@@ -132,6 +134,7 @@ namespace Opm
/// Capillary pressure.
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of nP capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of nP^2 derivative values,
/// array must be valid before calling.
@@ -140,6 +143,7 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void SaturationPropsFromDeck::capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const
{
@@ -147,12 +151,12 @@ namespace Opm
if (dpcds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i);
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalPc(s + np*i, pc + np*i);
funcForCell(cells[i]).evalPc(s + np*i, pc + np*i);
}
}
}
@@ -162,30 +166,84 @@ namespace Opm
/// Obtain the range of allowable saturation values.
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void SaturationPropsFromDeck::satRange(const int n,
const int* cells,
double* smin,
double* smax) const
{
const int np = phase_usage_.num_phases;
for (int i = 0; i < n; ++i) {
for (int p = 0; p < np; ++p) {
smin[np*i + p] = smin_[p];
smax[np*i + p] = smax_[p];
smin[np*i + p] = funcForCell(cells[i]).smin_[p];
smax[np*i + p] = funcForCell(cells[i]).smax_[p];
}
}
}
// Private methods below.
void SaturationPropsFromDeck::evalKr(const double* s, double* kr) const
// Map the cell number to the correct function set.
const SaturationPropsFromDeck::SatFuncSet&
SaturationPropsFromDeck::funcForCell(const int cell) const
{
if (phase_usage_.num_phases == 3) {
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// ----------- Methods of SatFuncSet below -----------
void SaturationPropsFromDeck::SatFuncSet::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg)
{
phase_usage = phase_usg;
const int samples = 200;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SaturationPropsFromDeck::SatFuncSet::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
@@ -203,18 +261,18 @@ namespace Opm
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage_.phase_used[Aqua]) {
int wpos = phase_usage_.phase_pos[Aqua];
int opos = phase_usage_.phase_pos[Liquid];
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
ASSERT(phase_usage_.phase_used[Vapour]);
int gpos = phase_usage_.phase_pos[Vapour];
int opos = phase_usage_.phase_pos[Liquid];
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
@@ -224,9 +282,9 @@ namespace Opm
}
void SaturationPropsFromDeck::evalKrDeriv(const double* s, double* kr, double* dkrds) const
void SaturationPropsFromDeck::SatFuncSet::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage_.num_phases;
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
@@ -256,9 +314,9 @@ namespace Opm
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage_.phase_used[Aqua]) {
int wpos = phase_usage_.phase_pos[Aqua];
int opos = phase_usage_.phase_pos[Liquid];
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
@@ -269,9 +327,9 @@ namespace Opm
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
ASSERT(phase_usage_.phase_used[Vapour]);
int gpos = phase_usage_.phase_pos[Vapour];
int opos = phase_usage_.phase_pos[Liquid];
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
@@ -286,20 +344,20 @@ namespace Opm
}
void SaturationPropsFromDeck::evalPc(const double* s, double* pc) const
void SaturationPropsFromDeck::SatFuncSet::evalPc(const double* s, double* pc) const
{
pc[phase_usage_.phase_pos[Liquid]] = 0.0;
if (phase_usage_.phase_used[Aqua]) {
int pos = phase_usage_.phase_pos[Aqua];
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage_.phase_used[Vapour]) {
int pos = phase_usage_.phase_pos[Vapour];
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SaturationPropsFromDeck::evalPcDeriv(const double* s, double* pc, double* dpcds) const
void SaturationPropsFromDeck::SatFuncSet::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
@@ -307,16 +365,16 @@ namespace Opm
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage_.num_phases;
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage_.phase_pos[Liquid]] = 0.0;
if (phase_usage_.phase_used[Aqua]) {
int pos = phase_usage_.phase_pos[Aqua];
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage_.phase_used[Vapour]) {
int pos = phase_usage_.phase_pos[Vapour];
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}

46
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@@ -23,6 +23,7 @@
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <vector>
namespace Opm
{
@@ -34,7 +35,9 @@ namespace Opm
SaturationPropsFromDeck();
/// Initialize from deck.
void init(const EclipseGridParser& deck);
/// global_cell maps from grid cells to their original logical Cartesian indices.
void init(const EclipseGridParser& deck,
const std::vector<int>& global_cell);
/// \return P, the number of phases.
int numPhases() const;
@@ -50,6 +53,7 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const;
@@ -64,6 +68,7 @@ namespace Opm
/// and is output in Fortran order (m_00 m_10 m_20 m01 ...)
void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const;
@@ -72,25 +77,36 @@ namespace Opm
/// \param[out] smin Array of nP minimum s values, array must be valid before calling.
/// \param[out] smax Array of nP maximum s values, array must be valid before calling.
void satRange(const int n,
const int* cells,
double* smin,
double* smax) const;
private:
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
private:
PhaseUsage phase_usage_;
utils::UniformTableLinear<double> krw_;
utils::UniformTableLinear<double> krow_;
utils::UniformTableLinear<double> pcow_;
utils::UniformTableLinear<double> krg_;
utils::UniformTableLinear<double> krog_;
utils::UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
class SatFuncSet
{
public:
void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
const SatFuncSet& funcForCell(const int cell) const;
};

4
opm/core/fluid/blackoil/SinglePvtDead.hpp
View File

@@ -74,8 +74,8 @@ namespace Opm
double* output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
utils::UniformTableLinear<double> one_over_B_;
utils::UniformTableLinear<double> viscosity_;
UniformTableLinear<double> one_over_B_;
UniformTableLinear<double> viscosity_;
};
}

200
opm/core/newwells.h
View File

@@ -20,69 +20,193 @@
#ifndef OPM_NEWWELLS_H_INCLUDED
#define OPM_NEWWELLS_H_INCLUDED
/**
* \file
*
* Main OPM-Core well data structure along with functions
* to create, populate and destroy it.
*/
#ifdef __cplusplus
extern "C" {
#endif
enum well_type { INJECTOR, PRODUCER };
enum control_type { BHP , RATE };
enum surface_component { WATER = 0, OIL = 1, GAS = 2 };
/** Well type indicates desired/expected well behaviour. */
enum WellType { INJECTOR, PRODUCER };
/** Type of well control equation or inequality constraint.
* BHP -> Well constrained by bottom-hole pressure target.
* RATE -> Well constrained by total reservoir volume flow rate.
*/
enum WellControlType { BHP , RATE };
/** Canonical component names and ordering. */
enum SurfaceComponent { WATER = 0, OIL = 1, GAS = 2 };
/* Control for single well */
struct WellControls {
int num;
int cpty;
enum control_type *type;
double *target;
int current;
/** Controls for a single well.
* Each control specifies a well rate or bottom-hole pressure. Only
* one control can be active at a time, indicated by current. The
* meaning of each control's target value depends on the control
* type, for BHP controls it is a pressure in Pascal, for RATE
* controls it is a volumetric rate in cubic(meter)/second. The
* active control as an equality constraint, whereas the
* non-active controls should be interpreted as inequality
* constraints (upper or lower bounds). For instance, a PRODUCER's BHP
* constraint defines a minimum acceptable bottom-hole pressure value
* for the well.
*/
struct WellControls
{
int num; /** Number of controls. */
enum WellControlType *type; /** Array of control types. */
double *target; /** Array of control targets. */
int current; /** Index of current active control. */
void *data; /** Internal management structure. */
};
struct Wells {
int number_of_wells;
int well_cpty;
int perf_cpty;
enum well_type *type;
double *depth_ref;
double *zfrac; /* Surface volume fraction
* (3*number_of_wells) */
int *well_connpos;
int *well_cells;
double *WI; /* Well index */
/** Data structure aggregating static information about all wells in a scenario. */
struct Wells
{
int number_of_wells; /** Number of wells. */
struct WellControls **ctrls; /* One struct for each well */
enum WellType *type; /** Array of well types. */
double *depth_ref; /** Array of well bhp reference depths. */
double *zfrac; /** Component volume fractions for each well, size is (3*number_of_wells).
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
*/
int *well_connpos; /** Array of indices into well_cells (and WI).
* For a well w, well_connpos[w] and well_connpos[w+1] yield
* start and one-beyond-end indices into the well_cells array
* for accessing w's perforation cell indices.
*/
int *well_cells; /** Array of perforation cell indices.
* Size is number of perforations (== well_connpos[number_of_wells]).
*/
double *WI; /** Well productivity index, same size and structure as well_cells. */
struct WellControls **ctrls; /** Well controls, one set of controls for each well. */
void *data; /** Internal management structure. */
};
struct CompletionData {
double *gpot; /* Gravity potential */
double *A; /* RB^{-1} */
double *phasemob; /* Phase mobility */
/** Data structure aggregating dynamic information about all wells in a scenario.
* All arrays in this structure contain data for each perforation,
* ordered the same as Wells::well_cells and Wells:WI. The array
* sizes are, respectively,
*
* gpot n*NP
* A n²*NP (matrix in column-major (i.e., Fortran) order).
* phasemob n*NP
*
* in which "n" denotes the number of active fluid phases (and
* constituent components) and "NP" is the total number of
* perforations, <CODE>well_connpos[ number_of_wells ]</CODE>.
*/
struct CompletionData
{
double *gpot; /** Gravity potentials. */
double *A; /** Volumes to surface-components matrix, A = RB^{-1}. */
double *phasemob; /** Phase mobilities. */
};
/**
* Construct a Wells object initially capable of managing a given
* number of wells and total number of well connections
* (perforations).
*
* Function add_well() is used to populate the Wells object. No
* reallocation occurs in function add_well() as long as the
* initially indicated capacities are sufficient. Call function
* destroy_wells() to dispose of the Wells object and its allocated
* memory resources.
*
* \param[in] nwells Expected number of wells in simulation scenario.
* Pass zero if the total number of wells is unknown.
*
* \param[in] nperf Expected total number of well connections
* (perforations) for all wells in simulation
* scenario. Pass zero if the total number of well
* connections is unknown.
*
* \return A valid Wells object with no wells if successful, and NULL
* otherwise.
*/
struct Wells *
wells_create(int nwells, int nperf);
create_wells(int nwells, int nperf);
/**
* Append a new well to an existing Wells object.
*
* Increments W->number_of_wells by one if successful. The new well
* does not include operational constraints. Such information is
* specified using function append_well_controls(). The current
* control index is set to -1 (invalid).
*
* \param[in] type Type of well.
* \param[in] depth_ref Reference depth for well's BHP.
* \param[in] nperf Number of perforations.
* \param[in] zfrac Injection fraction (three components) or NULL.
* \param[in] cells Grid cells in which well is perforated. Should
* ideally be track ordered.
* \param[in] WI Well production index per perforation, or NULL.
* \param[in,out] W Existing set of wells to which new well will
* be added.
*
* \return Non-zero (true) if successful and zero otherwise.
*/
int
wells_add(enum well_type type ,
double depth_ref,
int nperf ,
const double *zfrac , /* Injection fraction or NULL */
const int *cells ,
const double *WI , /* Well index per perf (or NULL) */
struct Wells *W );
add_well(enum WellType type ,
double depth_ref,
int nperf ,
const double *zfrac ,
const int *cells ,
const double *WI ,
struct Wells *W );
void
wells_destroy(struct Wells *W);
/**
* Append operational constraint to an existing well.
*
* Increments ctrl->num by one if successful. Introducing a new
* operational constraint does not affect the well's notion of the
* currently active constraint represented by ctrl->current.
*
* \param[in] type Control type.
* \param[in] target Target value for the control.
* \param[inout] ctrl Existing set of well controls.
* \return Non-zero (true) if successful and zero (false) otherwise.
*/
int
well_controls_append(enum control_type type ,
append_well_controls(enum WellControlType type ,
double target,
struct WellControls *ctrl );
/**
* Clear all controls from a well.
*
* Does not affect the control set capacity. */
void
well_controls_clear(struct WellControls *ctrl);
clear_well_controls(struct WellControls *ctrl);
/**
* Wells object destructor.
*
* Disposes of all resources managed by the Wells object.
*
* The Wells object must be built using function create_wells() and
* subsequently populated using function add_well().
*/
void
destroy_wells(struct Wells *W);
#ifdef __cplusplus
}

46
opm/core/pressure/tpfa/cfs_tpfa_residual.h
View File

@@ -54,15 +54,15 @@ struct cfs_tpfa_res_data {
struct cfs_tpfa_res_data *
cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
struct WellCompletions *wconn ,
int nphases);
cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
struct WellCompletions *wconn ,
int nphases);
void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
double dt,
struct cfs_tpfa_res_forces *forces,
const double *zc,
@@ -75,7 +75,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_forces *forces ,
int np ,
const double *trans ,
@@ -88,7 +88,7 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
double *wflux );
void
cfs_tpfa_res_fpress(struct UnstructuredGrid *G,
cfs_tpfa_res_fpress(struct UnstructuredGrid *G,
int np,
const double *htrans,
const double *pmobf,
@@ -100,19 +100,19 @@ cfs_tpfa_res_fpress(struct UnstructuredGrid *G,
#if 0
void
cfs_tpfa_retrieve_masstrans(struct UnstructuredGrid *G,
int np,
struct cfs_tpfa_data *h,
double *masstrans_f);
cfs_tpfa_retrieve_masstrans(struct UnstructuredGrid *G,
int np,
struct cfs_tpfa_data *h,
double *masstrans_f);
void
cfs_tpfa_retrieve_gravtrans(struct UnstructuredGrid *G,
int np,
struct cfs_tpfa_data *h,
double *gravtrans_f);
cfs_tpfa_retrieve_gravtrans(struct UnstructuredGrid *G,
int np,
struct cfs_tpfa_data *h,
double *gravtrans_f);
double
cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G,
cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G,
struct compr_quantities *cq,
const double *trans,
const double *porevol,
@@ -122,14 +122,14 @@ cfs_tpfa_impes_maxtime(struct UnstructuredGrid *G,
const double *gravity);
void
cfs_tpfa_expl_mass_transport(struct UnstructuredGrid *G,
well_t *W,
struct completion_data *wdata,
int np,
double dt,
const double *porevol,
struct cfs_tpfa_data *h,
double *surf_vol);
cfs_tpfa_expl_mass_transport(struct UnstructuredGrid *G,
well_t *W,
struct completion_data *wdata,
int np,
double dt,
const double *porevol,
struct cfs_tpfa_data *h,
double *surf_vol);
#endif

10
opm/core/transport/reorder/TransportModelTwophase.cpp
View File

@@ -549,7 +549,7 @@ namespace Opm
void TransportModelTwophase::solveGravityColumn(const std::vector<int>& cells)
int TransportModelTwophase::solveGravityColumn(const std::vector<int>& cells)
{
// Set up column gravflux.
const int nc = cells.size();
@@ -597,7 +597,7 @@ namespace Opm
THROW("In solveGravityColumn(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change);
}
// Repeat if necessary.
return num_iters + 1;
}
@@ -629,11 +629,13 @@ namespace Opm
saturation_ = &saturation[0];
// Solve on all columns.
int num_iters = 0;
for (std::vector<std::vector<int> >::size_type i = 0; i < columns.second.size(); i++) {
// std::cout << "==== new column" << std::endl;
solveGravityColumn(columns.second[i]);
num_iters += solveGravityColumn(columns.second[i]);
}
std::cout << "Gauss-Seidel column solver average iterations: "
<< double(num_iters)/double(columns.second.size()) << std::endl;
}
} // namespace Opm

2
opm/core/transport/reorder/TransportModelTwophase.hpp
View File

@@ -52,7 +52,7 @@ namespace Opm
void solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux);
void solveGravityColumn(const std::vector<int>& cells);
int solveGravityColumn(const std::vector<int>& cells);
void solveGravity(const std::pair<std::vector<int>, std::vector<std::vector<int> > >& columns,
const double* porevolume,
const double dt,

9
opm/core/utility/UniformTableLinear.hpp
View File

@@ -29,8 +29,6 @@
#include <opm/core/utility/ErrorMacros.hpp>
namespace Opm {
namespace utils {
/// @brief This class uses linear interpolation to compute the value
/// (and its derivative) of a function f sampled at uniform points.
@@ -253,7 +251,12 @@ namespace Opm {
return os;
}
} // namespace utils
namespace utils
{
using Opm::UniformTableLinear;
}
} // namespace Opm
#endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED

76
opm/core/utility/Units.hpp
View File

@@ -13,8 +13,8 @@
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
@@ -38,6 +38,7 @@
namespace Opm
{
namespace prefix
/// Conversion prefix for units.
{
const double micro = 1.0e-6;
const double milli = 1.0e-3;
@@ -49,68 +50,102 @@ namespace Opm
} // namespace prefix
namespace unit
/// Definition of various units.
/// All the units are defined in terms of international standard
/// units (SI). Example of use: We define a variable \c k which
/// gives a permeability. We want to set \c k to \f$1\,mD\f$.
/// \code
/// using namespace Opm::unit
/// double k = 0.001*darcy;
/// \endcode
/// We can also use one of the prefixes defined in Opm::prefix
/// \code
/// using namespace Opm::unit
/// using namespace Opm::prefix
/// double k = 1.0*milli*darcy;
/// \endcode
{
// Common powers
///\name Common powers
/// @{
inline double square(double v) { return v * v; }
inline double cubic (double v) { return v * v * v; }
/// @}
// --------------------------------------------------------------
// Basic (fundamental) units and conversions
// --------------------------------------------------------------
// Length:
/// \name Length
/// @{
const double meter = 1;
const double inch = 2.54 * prefix::centi*meter;
const double feet = 12 * inch;
/// @}
// Time:
/// \name Time
/// @{
const double second = 1;
const double minute = 60 * second;
const double hour = 60 * minute;
const double day = 24 * hour;
const double year = 365 * day;
/// @}
// Volume
const double stb = 0.158987294928 * cubic(meter);
/// \name Volume
/// @{
const double gallon = 231 * cubic(inch);
const double stb = 42 * gallon;
/// @}
// Mass:
/// \name Mass
/// @{
const double kilogram = 1;
// http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound
const double pound = 0.45359237 * kilogram;
/// @}
// --------------------------------------------------------------
// Standardised constants
// --------------------------------------------------------------
/// \name Standardised constant
/// @{
const double gravity = 9.80665 * meter/square(second);
/// @}
// --------------------------------------------------------------
// Derived units and conversions
// --------------------------------------------------------------
// Force:
/// \name Force
/// @{
const double Newton = kilogram*meter / square(second); // == 1
const double lbf = pound * gravity; // Pound-force
/// @}
// Pressure:
/// \name Pressure
/// @{
const double Pascal = Newton / square(meter); // == 1
const double barsa = 100000 * Pascal;
const double atm = 101325 * Pascal;
const double psia = lbf / square(inch);
/// @}
// Viscosity:
/// \name Viscosity
/// @{
const double Pas = Pascal * second; // == 1
const double Poise = prefix::deci*Pas;
/// @}
// Permeability:
//
// A porous medium with a permeability of 1 darcy permits a
// flow (flux) of 1 cm<63>/s of a fluid with viscosity 1 cP (1
// mPa<50>s) under a pressure gradient of 1 atm/cm acting across
// an area of 1 cm<63>.
//
/// \name Permeability
/// @{
///
/// A porous medium with a permeability of 1 darcy permits a
/// flow (flux) of \f$1\,cm^3/s\f$ of a fluid with viscosity
/// \f$1\,cP\f$ (\f$1\,mPa\cdot s\f$) under a pressure
/// gradient of \f$1\,atm/cm\f$ acting across an area of
/// \f$1\,cm^2\f$.
///
namespace perm_details {
const double p_grad = atm / (prefix::centi*meter);
const double area = square(prefix::centi*meter);
@@ -122,6 +157,7 @@ namespace Opm
// == 9.869232667160130e-13 [m^2]
}
const double darcy = perm_details::darcy;
/// @}
// Unit conversion support.
//

34
opm/core/utility/buildUniformMonotoneTable.hpp
View File

@@ -24,27 +24,25 @@
#include <opm/core/utility/UniformTableLinear.hpp>
namespace Opm {
namespace utils {
template <typename T>
void buildUniformMonotoneTable(const std::vector<double>& xv,
const std::vector<T>& yv,
const int samples,
UniformTableLinear<T>& table)
{
MonotCubicInterpolator interp(xv, yv);
std::vector<T> uniform_yv(samples);
double xmin = xv[0];
double xmax = xv.back();
for (int i = 0; i < samples; ++i) {
double w = double(i)/double(samples - 1);
double x = (1.0 - w)*xmin + w*xmax;
uniform_yv[i] = interp(x);
}
table = UniformTableLinear<T>(xmin, xmax, uniform_yv);
template <typename T>
void buildUniformMonotoneTable(const std::vector<double>& xv,
const std::vector<T>& yv,
const int samples,
UniformTableLinear<T>& table)
{
MonotCubicInterpolator interp(xv, yv);
std::vector<T> uniform_yv(samples);
double xmin = xv[0];
double xmax = xv.back();
for (int i = 0; i < samples; ++i) {
double w = double(i)/double(samples - 1);
double x = (1.0 - w)*xmin + w*xmax;
uniform_yv[i] = interp(x);
}
table = UniformTableLinear<T>(xmin, xmax, uniform_yv);
}
} // namespace utils
} // namespace Opm

163
opm/core/utility/newwells.c
View File

@@ -41,14 +41,70 @@
#include <stdlib.h>
#include <string.h>
struct WellControlMgmt {
int cpty;
};
struct WellMgmt {
int well_cpty;
int perf_cpty;
};
static void
destroy_ctrl_mgmt(struct WellControlMgmt *m)
{
free(m);
}
static struct WellControlMgmt *
create_ctrl_mgmt(void)
{
struct WellControlMgmt *m;
m = malloc(1 * sizeof *m);
if (m != NULL) {
m->cpty = 0;
}
return m;
}
static void
destroy_well_mgmt(struct WellMgmt *m)
{
free(m);
}
static struct WellMgmt *
create_well_mgmt(void)
{
struct WellMgmt *m;
m = malloc(1 * sizeof *m);
if (m != NULL) {
m->well_cpty = 0;
m->perf_cpty = 0;
}
return m;
}
/* ---------------------------------------------------------------------- */
static void
well_controls_destroy(struct WellControls *ctrl)
/* ---------------------------------------------------------------------- */
{
if (ctrl != NULL) {
free(ctrl->target);
free(ctrl->type);
destroy_ctrl_mgmt(ctrl->data);
free (ctrl->target);
free (ctrl->type);
}
free(ctrl);
@@ -67,10 +123,16 @@ well_controls_create(void)
if (ctrl != NULL) {
/* Initialise empty control set */
ctrl->num = 0;
ctrl->cpty = 0;
ctrl->type = NULL;
ctrl->target = NULL;
ctrl->current = -1;
ctrl->data = create_ctrl_mgmt();
if (ctrl->data == NULL) {
well_controls_destroy(ctrl);
ctrl = NULL;
}
}
return ctrl;
@@ -85,6 +147,8 @@ well_controls_reserve(int nctrl, struct WellControls *ctrl)
int c, ok;
void *type, *target;
struct WellControlMgmt *m;
type = realloc(ctrl->type , nctrl * sizeof *ctrl->type );
target = realloc(ctrl->target, nctrl * sizeof *ctrl->target);
@@ -93,12 +157,13 @@ well_controls_reserve(int nctrl, struct WellControls *ctrl)
if (target != NULL) { ctrl->target = target; ok++; }
if (ok == 2) {
for (c = ctrl->cpty; c < nctrl; c++) {
m = ctrl->data;
for (c = m->cpty; c < nctrl; c++) {
ctrl->type [c] = BHP;
ctrl->target[c] = -1.0;
}
ctrl->cpty = nctrl;
m->cpty = nctrl;
}
return ok == 2;
@@ -160,7 +225,11 @@ initialise_new_wells(int nwells, struct Wells *W)
{
int ok, w;
for (w = W->well_cpty; w < nwells; w++) {
struct WellMgmt *m;
m = W->data;
for (w = m->well_cpty; w < nwells; w++) {
W->type [w] = PRODUCER;
W->depth_ref[w] = -1.0;
@@ -171,7 +240,7 @@ initialise_new_wells(int nwells, struct Wells *W)
W->well_connpos[w + 1] = W->well_connpos[w];
}
for (w = W->well_cpty, ok = 1; ok && (w < nwells); w++) {
for (w = m->well_cpty, ok = 1; ok && (w < nwells); w++) {
W->ctrls[w] = well_controls_create();
ok = W->ctrls[w] != NULL;
@@ -194,7 +263,11 @@ initialise_new_perfs(int nperf, struct Wells *W)
{
int k;
for (k = W->perf_cpty; k < nperf; k++) {
struct WellMgmt *m;
m = W->data;
for (k = m->perf_cpty; k < nperf; k++) {
W->well_cells[k] = -1 ;
W->WI [k] = 0.0;
}
@@ -208,12 +281,16 @@ wells_reserve(int nwells, int nperf, struct Wells *W)
{
int ok;
assert (nwells >= W->well_cpty);
assert (nperf >= W->perf_cpty);
struct WellMgmt *m;
m = W->data;
assert (nwells >= m->well_cpty);
assert (nperf >= m->perf_cpty);
ok = 1;
if (nwells > W->well_cpty) {
if (nwells > m->well_cpty) {
ok = wells_allocate(nwells, W);
if (ok) {
@@ -221,16 +298,16 @@ wells_reserve(int nwells, int nperf, struct Wells *W)
}
if (ok) {
W->well_cpty = nwells;
m->well_cpty = nwells;
}
}
if (ok && (nperf > W->perf_cpty)) {
if (ok && (nperf > m->perf_cpty)) {
ok = perfs_allocate(nperf, W);
if (ok) {
initialise_new_perfs(nperf, W);
W->perf_cpty = nperf;
m->perf_cpty = nperf;
}
}
@@ -245,7 +322,7 @@ wells_reserve(int nwells, int nperf, struct Wells *W)
/* ---------------------------------------------------------------------- */
struct Wells *
wells_create(int nwells, int nperf)
create_wells(int nwells, int nperf)
/* ---------------------------------------------------------------------- */
{
int ok;
@@ -255,8 +332,6 @@ wells_create(int nwells, int nperf)
if (W != NULL) {
W->number_of_wells = 0;
W->well_cpty = 0;
W->perf_cpty = 0;
W->type = NULL;
W->depth_ref = NULL;
@@ -268,7 +343,9 @@ wells_create(int nwells, int nperf)
W->ctrls = NULL;
ok = W->well_connpos != NULL;
W->data = create_well_mgmt();
ok = (W->well_connpos != NULL) && (W->data != NULL);
if (ok) {
W->well_connpos[0] = 0;
@@ -278,7 +355,7 @@ wells_create(int nwells, int nperf)
}
if (! ok) {
wells_destroy(W);
destroy_wells(W);
W = NULL;
}
}
@@ -289,16 +366,22 @@ wells_create(int nwells, int nperf)
/* ---------------------------------------------------------------------- */
void
wells_destroy(struct Wells *W)
destroy_wells(struct Wells *W)
/* ---------------------------------------------------------------------- */
{
int w;
struct WellMgmt *m;
if (W != NULL) {
for (w = 0; w < W->well_cpty; w++) {
m = W->data;
for (w = 0; w < m->well_cpty; w++) {
well_controls_destroy(W->ctrls[w]);
}
destroy_well_mgmt(m);
free(W->ctrls);
free(W->WI);
free(W->well_cells);
@@ -331,26 +414,30 @@ alloc_size(int n, int a, int cpty)
/* ---------------------------------------------------------------------- */
int
wells_add(enum well_type type ,
double depth_ref,
int nperf ,
const double *zfrac , /* Injection fraction or NULL */
const int *cells ,
const double *WI , /* Well index per perf (or NULL) */
struct Wells *W )
add_well(enum WellType type ,
double depth_ref,
int nperf ,
const double *zfrac , /* Injection fraction or NULL */
const int *cells ,
const double *WI , /* Well index per perf (or NULL) */
struct Wells *W )
/* ---------------------------------------------------------------------- */
{
int ok, nw, nperf_tot, off;
int nwalloc, nperfalloc;
struct WellMgmt *m;
nw = W->number_of_wells;
nperf_tot = W->well_connpos[nw];
ok = (nw < W->well_cpty) && (nperf_tot + nperf <= W->perf_cpty);
m = W->data;
ok = (nw < m->well_cpty) && (nperf_tot + nperf <= m->perf_cpty);
if (! ok) {
nwalloc = alloc_size(nw , 1 , W->well_cpty);
nperfalloc = alloc_size(nperf_tot, nperf, W->perf_cpty);
nwalloc = alloc_size(nw , 1 , m->well_cpty);
nperfalloc = alloc_size(nperf_tot, nperf, m->perf_cpty);
ok = wells_reserve(nwalloc, nperfalloc, W);
}
@@ -386,19 +473,22 @@ wells_add(enum well_type type ,
/* ---------------------------------------------------------------------- */
int
well_controls_append(enum control_type type ,
append_well_controls(enum WellControlType type ,
double target,
struct WellControls *ctrl )
/* ---------------------------------------------------------------------- */
{
int ok, alloc;
struct WellControlMgmt *m;
assert (ctrl != NULL);
ok = ctrl->num < ctrl->cpty;
m = ctrl->data;
ok = ctrl->num < m->cpty;
if (! ok) {
alloc = alloc_size(ctrl->num, 1, ctrl->cpty);
alloc = alloc_size(ctrl->num, 1, m->cpty);
ok = well_controls_reserve(alloc, ctrl);
}
@@ -407,9 +497,6 @@ well_controls_append(enum control_type type ,
ctrl->target[ctrl->num] = target;
ctrl->num += 1;
/* TODO: Review this: */
ctrl->current = 0;
}
return ok;
@@ -418,7 +505,7 @@ well_controls_append(enum control_type type ,
/* ---------------------------------------------------------------------- */
void
well_controls_clear(struct WellControls *ctrl)
clear_well_controls(struct WellControls *ctrl)
/* ---------------------------------------------------------------------- */
{
if (ctrl != NULL) {

5
tutorials/Makefile.am
View File

@@ -13,3 +13,8 @@ if UMFPACK
noinst_PROGRAMS += tutorial2
tutorial2_SOURCES = tutorial2.cpp
endif
if UMFPACK
noinst_PROGRAMS += tutorial3
tutorial3_SOURCES = tutorial3.cpp
endif

36
tutorials/tutorial1.cpp
View File

@@ -26,43 +26,41 @@
#include "config.h"
#endif // HAVE_CONFIG_H
/// \page tutorial1 A simple carthesian grid
/// This tutorial explains how to construct a simple carthesian grid.\n\n
/// We construct a 2x2 two dimensional carthesian grid with 4 blocks of equal size.
/// \page tutorial1 A simple cartesian grid
/// This tutorial explains how to construct a simple cartesian grid.
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <cassert>
#include <cstddef>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
// ----------------- Main program -----------------
/// \page tutorial1
/// \page tutorial1
/// \section commentedsource1 Commented source code:
/// \code
int main()
{
/// \endcode
/// \page tutorial1
/// By setting <code>nz = 1</code>, we make the grid two dimensional
/// We set the number of blocks in each direction.
/// \code
int nx = 3;
int nx = 4;
int ny = 3;
int nz = 1;
int nz = 2;
/// \endcode
/// The size of each block is 1x1x1. We use standard units (SI)
/// The size of each block is 1x1x1. The default units are allways the
/// standard units (SI). But other units can easily be dealt with, see Opm::unit.
/// \code
double dx = 1.0;
double dy = 1.0;
double dz = 1.0;
/// \endcode
/// \page tutorial1
/// One of the constructors of the class Opm::GridManager takes <code>nx,ny,nz,dx,dy,dz</code>
/// and construct the corresponding carthesian grid.
/// \endcode
/// \page tutorial1
/// One of the constructors of the class Opm::GridManager takes <code>nx,ny,nz,dx,dy,dz</code>
/// and construct the corresponding cartesian grid.
/// \code
Opm::GridManager grid(nx, ny, nz, dx, dy, dz);
/// \endcode
@@ -84,10 +82,10 @@ int main()
}
/// \endcode
/// \page tutorial1
/// We read the the vtu output file in \a Paraview and obtain the following grid.
/// We read the vtu output file in \a Paraview and obtain the following grid.
/// \image html tutorial1.png
/// \page tutorial1
/// \section sourcecode Source code.
/// \include tutorial1.cpp
/// \page tutorial1
/// \section completecode1 Complete source code:
/// \include tutorial1.cpp

53
tutorials/tutorial2.cpp
View File

@@ -18,14 +18,14 @@
*/
/// \page tutorial2 Flow Solver for a single phase
/// \page tutorial2 Flow Solver for a single phase
/// \details The flow equations consist of the mass conservation equation
/// \f[\nabla\cdot u=q\f] and the Darcy law \f[u=-\frac{1}{\mu}K\nabla p.\f] Here,
/// \f$u\f$ denotes the velocity and \f$p\f$ the pressure. The permeability tensor is
/// given by \f$K\f$ and \f$\mu\f$ denotes the viscosity.
///
/// We solve the flow equations for a carthesian grid and we set the source term
/// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal
/// given by \f$K\f$ and \f$\mu\f$ denotes the viscosity.
///
/// We solve the flow equations for a cartesian grid and we set the source term
/// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal
/// with opposite sign (inflow equal to outflow).
@@ -36,25 +36,23 @@
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <cassert>
#include <cstddef>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
/// \page tutorial2
/// \section commentedcode Commented code:
/// \section commentedcode2 Commented code:
/// \code
int main()
{
/// \endcode
/// \page tutorial2
/// We construct a carthesian grid
/// We construct a cartesian grid
/// \code
int dim = 3;
int nx = 40;
@@ -71,19 +69,25 @@ int main()
int num_faces = grid.c_grid()->number_of_faces;
/// \endcode
/// \page tutorial2
/// We define the viscosity (unit: cP).
/// \details
/// We define a fluid viscosity equal to \f$1\,cP\f$. We use
/// the namespaces Opm::unit
/// and Opm::prefix to deal with the units.
/// \code
double mu = 1.0;
using namespace Opm::unit;
using namespace Opm::prefix;
double mu = 1.0*centi*Poise;
/// \endcode
/// \page tutorial2
/// We define the permeability (unit: mD).
/// \details
/// We define a permeability equal to \f$100\,mD\f$.
/// \code
double k = 100.0;
double k = 100.0*milli*darcy;
/// \endcode
/// \page tutorial2
/// \details
/// We set up the permeability tensor and compute the mobility for each cell.
/// The permeability tensor is flattened in a vector.
/// The resulting permeability matrix is flattened in a vector.
/// \code
std::vector<double> permeability(num_cells*dim*dim, 0.);
std::vector<double> mob(num_cells);
@@ -96,7 +100,8 @@ int main()
/// \endcode
/// \page tutorial2
/// We choose the UMFPACK linear solver for the pressure solver.
/// \details
/// We take UMFPACK as the linear solver for the pressure solver (This library has therefore to be installed.)
/// \code
Opm::LinearSolverUmfpack linsolver;
/// \endcode
@@ -115,7 +120,7 @@ int main()
src[num_cells-1] = -100.;
/// \endcode
/// \page tutorial2
/// \details We set up the boundary conditions. We do not modify them.
/// \details We set up the boundary conditions. We do not modify them.
/// By default, we obtain no outflow boundary conditions.
/// \code
/// \code
@@ -133,21 +138,21 @@ int main()
std::vector<double> faceflux(num_faces);
std::vector<double> well_bhp;
std::vector<double> well_flux;
/// \endcode
/// \endcode
/// \page tutorial2
/// \details
/// We declare the gravity term which is required by the pressure solver (see
/// Opm::IncompTpfa.solve()). In the absence of gravity, an empty vector is required.
/// \code
std::vector<double> omega;
std::vector<double> omega;
/// \endcode
/// \page tutorial2
/// \details
/// We declare the wdp term which is required by the pressure solver (see
/// Opm::IncompTpfa.solve()). In the absence of wells, an empty vector is required.
/// \code
std::vector<double> wdp;
std::vector<double> wdp;
/// \endcode
/// \page tutorial2
@@ -170,10 +175,10 @@ int main()
}
/// \endcode
/// \page tutorial2
/// We read the the vtu output file in \a Paraview and obtain the following pressure
/// We read the vtu output file in \a Paraview and obtain the following pressure
/// distribution. \image html tutorial2.png
/// \page tutorial2
/// \section sourcecode Complete source code.
/// \include tutorial2.cpp
/// \section completecode2 Complete source code:
/// \include tutorial2.cpp

521
tutorials/tutorial3.cpp Normal file
View File

@@ -0,0 +1,521 @@
/*
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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cassert>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/transport/transport_source.h>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/NormSupport.hpp>
#include <opm/core/transport/ImplicitAssembly.hpp>
#include <opm/core/transport/ImplicitTransport.hpp>
#include <opm/core/transport/JacobianSystem.hpp>
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/TwophaseState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
/// \page tutorial3 Multiphase flow
/// The Darcy law gives
/// \f[u_\alpha= -\frac1{\mu_\alpha} K_\alpha\nabla p_\alpha\f]
/// where \f$\mu_\alpha\f$ and \f$K_\alpha\f$ represent the viscosity
/// and the permeability tensor for each phase \f$\alpha\f$. In the two phase
/// case, we have either \f$\alpha=w\f$ or \f$\alpha=o\f$.
/// In this tutorial, we do not take into account capillary pressure so that
/// \f$p=p_w=p_o\f$ and gravity
/// effects. We denote by \f$K\f$ the absolute permeability tensor and each phase
/// permeability is defined through its relative permeability by the expression
/// \f[K_\alpha=k_{r\alpha}K.\f]
/// The phase mobility are defined as
/// \f[\lambda_\alpha=\frac{k_{r\alpha}}{\mu_\alpha}\f]
/// so that the Darcy law may be rewritten as
/// \f[u_\alpha= -\lambda_\alpha K\nabla p.\f]
/// The conservation of mass for each phase writes:
/// \f[\frac{\partial}{\partial t}(\phi\rho_\alpha s_\alpha)+\nabla\cdot(\rho_\alpha u_\alpha)=q_\alpha\f]
/// where \f$s_\alpha\f$ denotes the saturation of the phase \f$\alpha\f$ and \f$q_\alpha\f$ is a source term. Let
/// us consider a two phase flow with oil and water. We assume that the phases are incompressible. Since
/// \f$s_w+s_o=1\f$, we get
/// \f[ \nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o}.\f]
/// where we define
/// \f[u=u_w+u_o.\f]
/// Let the total mobility be equal to
/// \f[\lambda=\lambda_w+\lambda_o\f]
/// Then, we have
/// \f[u=-\lambda K\nabla p.\f]
/// The set of equations
/// \f[\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o},\quad u=-\lambda K\nabla p.\f]
/// is referred to as the <strong>pressure equation</strong>. We introduce
/// the fractional flow \f$f_w\f$
/// as
/// \f[f_w=\frac{\lambda_w}{\lambda_w+\lambda_o}\f]
/// and obtain
/// \f[\phi\frac{\partial s}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}\f]
/// which is referred to as the <strong>transport equation</strong>. The pressure and
/// transport equation are coupled. In this tutorial, we implement a splitting scheme,
/// where, at each time step, we decouple the two equations. We solve first
/// the pressure equation and then update the water saturation by solving
/// the transport equation assuming that \f$u\f$ is constant in time in the time step
/// interval we are considering.
/// \page tutorial3
/// \section commentedcode3 Commented code:
/// \page tutorial3
/// \details
/// We define a class which computes mobility, capillary pressure and
/// the minimum and maximum saturation value for each cell.
/// \code
class TwophaseFluid
{
public:
/// \endcode
/// \page tutorial3
/// \details Constructor operator. Takes in the fluid properties defined
/// \c props
/// \code
TwophaseFluid(const Opm::IncompPropertiesInterface& props);
/// \endcode
/// \page tutorial3
/// \details Density for each phase.
/// \code
double density(int phase) const;
/// \endcode
/// \page tutorial3
/// \details Computes the mobility and its derivative with respect to saturation
/// for a given cell \c c and saturation \c sat. The template classes \c Sat,
/// \c Mob, \c DMob are typically arrays. By using templates, we do not have to
/// investigate how these array objects are implemented
/// (as long as they have an \c operator[] method).
/// \code
template <class Sat, class Mob, class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const;
/// \endcode
/// \page tutorial3
/// \details Computes the capillary pressure and its derivative with respect
/// to saturation
/// for a given cell \c c and saturation \c sat.
/// \code
template <class Sat, class Pcap, class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const;
/// \endcode
/// \page tutorial3
/// \details Returns the minimum permitted saturation.
/// \code
double s_min(int c) const;
/// \endcode
/// \page tutorial3
/// \details Returns the maximum permitted saturation
/// \code
double s_max(int c) const;
/// \endcode
/// \page tutorial3
/// \details Private variables
/// \code
private:
const Opm::IncompPropertiesInterface& props_;
std::vector<double> smin_;
std::vector<double> smax_;
};
/// \endcode
/// \page tutorial3
/// \details We set up the transport model.
/// \code
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
/// \endcode
/// \page tutorial3
/// \details
/// The transport equation is nonlinear. We use an implicit transport solver
/// which implements a Newton-Raphson solver.
/// We define the format of the objects
/// which will be used by the solver.
/// \code
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return AccumulationNorm <Vector, MaxAbs>::norm(v);
}
};
/// \endcode
/// \page tutorial3
/// \details
/// We set up the solver.
/// \code
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
/// \endcode
/// \page tutorial3
/// \details
/// Main function
/// \code
int main ()
{
/// \endcode
/// \page tutorial3
/// \details
/// We define the grid. A cartesian grid with 1200 cells.
/// \code
int dim = 3;
int nx = 20;
int ny = 20;
int nz = 1;
double dx = 10.;
double dy = 10.;
double dz = 10.;
using namespace Opm;
GridManager grid(nx, ny, nz, dx, dy, dz);
int num_cells = grid.c_grid()->number_of_cells;
/// \endcode
/// \page tutorial3
/// \details
/// We define the properties of the fluid.\n
/// Number of phases.
/// \code
int num_phases = 2;
using namespace unit;
using namespace prefix;
/// \endcode
/// \page tutorial3
/// \details density vector (one component per phase).
/// \code
std::vector<double> rho(2, 1000.);
/// \endcode
/// \page tutorial3
/// \details viscosity vector (one component per phase).
/// \code
std::vector<double> mu(2, 1.*centi*Poise);
/// \endcode
/// \page tutorial3
/// \details porosity and permeability of the rock.
/// \code
double porosity = 0.5;
double k = 10*milli*darcy;
/// \endcode
/// \page tutorial3
/// \details We define the relative permeability function. We use a basic fluid
/// description and set this function to be linear. For more realistic fluid, the
/// saturation function is given by the data.
/// \code
SaturationPropsBasic::RelPermFunc rel_perm_func;
rel_perm_func = SaturationPropsBasic::Linear;
/// \endcode
/// \page tutorial3
/// \details We construct a basic fluid with the properties we have defined above.
/// Each property is constant and hold for all cells.
/// \code
IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu,
porosity, k, dim, num_cells);
TwophaseFluid fluid(props);
/// \endcode
/// \page tutorial3
/// \details Gravity parameters. Here, we set zero gravity.
/// \code
const double *grav = 0;
std::vector<double> omega;
/// \endcode
/// \page tutorial3
/// \details We set up the pressure solver.
/// \code
LinearSolverUmfpack linsolver;
IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver);
/// \endcode
/// \page tutorial3
/// \details We set up the source term
/// \code
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \endcode
/// \page tutorial3
/// \details We set up the wells. Here, there are no well and we let them empty.
/// \code
std::vector<double> empty_wdp;
std::vector<double> empty_well_bhp;
std::vector<double> empty_well_flux;
/// \endcode
/// \page tutorial3
/// \details We set up the source term for the transport solver.
/// \code
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);
}
}
/// \endcode
/// \page tutorial3
/// \details We compute the pore volume
/// \code
std::vector<double> porevol;
computePorevolume(*grid.c_grid(), props, porevol);
/// \endcode
/// \page tutorial3
/// \details We set up the transport solver.
/// \code
const bool guess_old_solution = true;
TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution);
TransportSolver tsolver(model);
/// \endcode
/// \page tutorial3
/// \details Time integration parameters
/// \code
double dt = 0.1*day;
int num_time_steps = 20;
/// \endcode
/// \page tutorial3
/// \details Control paramaters for the implicit solver.
/// \code
ImplicitTransportDetails::NRReport rpt;
ImplicitTransportDetails::NRControl ctrl;
/// \endcode
/// \page tutorial3
/// \details We define a vector which contains all cell indexes. We use this
/// vector to set up parameters on the whole domains.
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
/// \page tutorial3
/// \details We set up the boundary conditions. Letting bcs empty is equivalent
/// to no flow boundary conditions.
/// \code
FlowBCManager bcs;
/// \endcode
/// \page tutorial3
/// \details
/// Linear solver init.
/// \code
using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
/// \endcode
/// \page tutorial3
/// \details
/// We initialise water saturation to minimum everywhere.
/// \code
TwophaseState state;
state.init(*grid.c_grid());
state.setWaterSat(allcells, props, TwophaseState::MinSat);
/// \endcode
/// \page tutorial3
/// \details We introduce a vector which contains the total mobility
/// on all cells.
/// \code
std::vector<double> totmob;
/// \endcode
/// \page tutorial3
/// \details This string will contain the name of a VTK output vector.
/// \code
std::ostringstream vtkfilename;
/// \endcode
/// \page tutorial3
/// \details Loop over the time steps.
/// \code
for (int i = 0; i < num_time_steps; ++i) {
/// \endcode
/// \page tutorial3
/// \details Compute the total mobility. It is needed by the pressure solver
/// \code
computeTotalMobility(props, allcells, state.saturation(), totmob);
/// \endcode
/// \page tutorial3
/// \details Solve the pressure equation
/// \code
psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(),
state.pressure(), state.faceflux(), empty_well_bhp,
empty_well_flux);
/// \endcode
/// \page tutorial3
/// \details Transport solver
/// \code
tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt);
/// \endcode
/// \page tutorial3
/// \details Write the output to file.
/// \code
vtkfilename.str("");
vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
}
}
/// \endcode
/// \page tutorial3
/// \details Implementation of the TwophaseFluid class.
/// \code
TwophaseFluid::TwophaseFluid(const Opm::IncompPropertiesInterface& props)
: props_(props),
smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases())
{
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 TwophaseFluid::density(int phase) const
{
return props_.density()[phase];
}
template <class Sat,
class Mob,
class DMob>
void TwophaseFluid::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];
/// \endcode
/// \page tutorial3
/// \details We
/// 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.
/// \code
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>
void TwophaseFluid::pc(int /*c */, const Sat& /* s*/, Pcap& pcap, DPcap& dpcap) const
{
pcap = 0.;
dpcap = 0.;
}
double TwophaseFluid::s_min(int c) const
{
return smin_[2*c + 0];
}
double TwophaseFluid::s_max(int c) const
{
return smax_[2*c + 0];
}
/// \endcode
/// \page tutorial3
/// \section results3 Results.
/// <TABLE>
/// <TR>
/// <TD> \image html tutorial3-000.png </TD>
/// <TD> \image html tutorial3-005.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial3-010.png </TD>
/// <TD> \image html tutorial3-015.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial3-019.png </TD>
/// <TD> </TD>
/// </TR>
/// </TABLE>
/// \page tutorial3
/// \details
/// \section completecode3 Complete source code:
/// \include tutorial3.cpp
/// \include generate_doc_figures.py