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

Merge pull request #730 from atgeirr/compute-tof

Browse Source 
More features for flow diagnostics
This commit is contained in:
Bård Skaflestad 2015-01-23 10:40:43 +01:00
commit 114b64f23d
9 changed files with 1147 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists_files.cmake
View File

@ -55,6 +55,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/pressure/CompressibleTpfa.cpp opm/core/pressure/CompressibleTpfa.cpp
opm/core/pressure/FlowBCManager.cpp opm/core/pressure/FlowBCManager.cpp
opm/core/pressure/IncompTpfa.cpp opm/core/pressure/IncompTpfa.cpp
opm/core/pressure/IncompTpfaSinglePhase.cpp
opm/core/pressure/cfsh.c opm/core/pressure/cfsh.c
opm/core/pressure/flow_bc.c opm/core/pressure/flow_bc.c
opm/core/pressure/fsh.c opm/core/pressure/fsh.c
@ -83,6 +84,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/props/BlackoilPropertiesFromDeck.cpp opm/core/props/BlackoilPropertiesFromDeck.cpp
opm/core/props/IncompPropertiesBasic.cpp opm/core/props/IncompPropertiesBasic.cpp
opm/core/props/IncompPropertiesFromDeck.cpp opm/core/props/IncompPropertiesFromDeck.cpp
opm/core/props/IncompPropertiesSinglePhase.cpp
opm/core/props/pvt/BlackoilPvtProperties.cpp opm/core/props/pvt/BlackoilPvtProperties.cpp
opm/core/props/pvt/PvtPropertiesBasic.cpp opm/core/props/pvt/PvtPropertiesBasic.cpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
@ -110,6 +112,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/simulator/SimulatorTimer.cpp opm/core/simulator/SimulatorTimer.cpp
opm/core/tof/AnisotropicEikonal.cpp opm/core/tof/AnisotropicEikonal.cpp
opm/core/tof/DGBasis.cpp opm/core/tof/DGBasis.cpp
opm/core/tof/FlowDiagnostics.cpp
opm/core/tof/TofReorder.cpp opm/core/tof/TofReorder.cpp
opm/core/tof/TofDiscGalReorder.cpp opm/core/tof/TofDiscGalReorder.cpp
opm/core/transport/TransportSolverTwophaseInterface.cpp opm/core/transport/TransportSolverTwophaseInterface.cpp
@ -162,6 +165,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_ug.cpp tests/test_ug.cpp
tests/test_cubic.cpp tests/test_cubic.cpp
tests/test_event.cpp tests/test_event.cpp
tests/test_flowdiagnostics.cpp
tests/test_nonuniformtablelinear.cpp tests/test_nonuniformtablelinear.cpp
tests/test_sparsevector.cpp tests/test_sparsevector.cpp
tests/test_sparsetable.cpp tests/test_sparsetable.cpp
@ -302,6 +306,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/pressure/CompressibleTpfa.hpp opm/core/pressure/CompressibleTpfa.hpp
opm/core/pressure/FlowBCManager.hpp opm/core/pressure/FlowBCManager.hpp
opm/core/pressure/IncompTpfa.hpp opm/core/pressure/IncompTpfa.hpp
opm/core/pressure/IncompTpfaSinglePhase.hpp
opm/core/pressure/flow_bc.h opm/core/pressure/flow_bc.h
opm/core/pressure/fsh.h opm/core/pressure/fsh.h
opm/core/pressure/fsh_common_impl.h opm/core/pressure/fsh_common_impl.h
@ -335,6 +340,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/IncompPropertiesInterface.hpp opm/core/props/IncompPropertiesInterface.hpp
opm/core/props/IncompPropertiesShadow.hpp opm/core/props/IncompPropertiesShadow.hpp
opm/core/props/IncompPropertiesShadow_impl.hpp opm/core/props/IncompPropertiesShadow_impl.hpp
opm/core/props/IncompPropertiesSinglePhase.hpp
opm/core/props/phaseUsageFromDeck.hpp opm/core/props/phaseUsageFromDeck.hpp
opm/core/props/pvt/BlackoilPvtProperties.hpp opm/core/props/pvt/BlackoilPvtProperties.hpp
opm/core/props/pvt/PvtPropertiesBasic.hpp opm/core/props/pvt/PvtPropertiesBasic.hpp
@ -379,6 +385,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/initStateEquil_impl.hpp opm/core/simulator/initStateEquil_impl.hpp
opm/core/tof/AnisotropicEikonal.hpp opm/core/tof/AnisotropicEikonal.hpp
opm/core/tof/DGBasis.hpp opm/core/tof/DGBasis.hpp
opm/core/tof/FlowDiagnostics.hpp
opm/core/tof/TofReorder.hpp opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp opm/core/tof/TofDiscGalReorder.hpp
opm/core/transport/TransportSolverTwophaseInterface.hpp opm/core/transport/TransportSolverTwophaseInterface.hpp

241
examples/compute_tof.cpp
View File

@ -23,8 +23,6 @@
#include "config.h" #include "config.h"
#endif // HAVE_CONFIG_H #endif // HAVE_CONFIG_H
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp> #include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h> #include <opm/core/wells.h>
@ -32,18 +30,15 @@
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp> #include <opm/core/utility/SparseTable.hpp>
#include <opm/core/utility/StopWatch.hpp> #include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp> #include <opm/core/props/IncompPropertiesSinglePhase.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp> #include <opm/core/pressure/IncompTpfaSinglePhase.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/tof/TofReorder.hpp> #include <opm/core/tof/TofReorder.hpp>
#include <opm/core/tof/TofDiscGalReorder.hpp> #include <opm/core/tof/TofDiscGalReorder.hpp>
@ -70,23 +65,88 @@ namespace
} }
} }
void buildTracerheadsFromWells(const Wells* wells, void buildTracerheadsFromWells(const Wells& wells,
const bool trace_injectors,
Opm::SparseTable<int>& tracerheads) Opm::SparseTable<int>& tracerheads)
{ {
if (wells == 0) {
return;
}
tracerheads.clear(); tracerheads.clear();
const int num_wells = wells->number_of_wells; const int num_wells = wells.number_of_wells;
const WellType wanted_type = trace_injectors ? INJECTOR : PRODUCER;
for (int w = 0; w < num_wells; ++w) { for (int w = 0; w < num_wells; ++w) {
if (wells->type[w] != INJECTOR) { if (wells.type[w] != wanted_type) {
continue; continue;
} }
tracerheads.appendRow(wells->well_cells + wells->well_connpos[w], tracerheads.appendRow(wells.well_cells + wells.well_connpos[w],
wells->well_cells + wells->well_connpos[w + 1]); wells.well_cells + wells.well_connpos[w + 1]);
} }
} }
void setBhpWells(Wells& wells)
{
const int num_wells = wells.number_of_wells;
for (int w = 0; w < num_wells; ++w) {
WellControls* ctrl = wells.ctrls[w];
const double target = (wells.type[w] == INJECTOR) ? 200*Opm::unit::barsa : 100*Opm::unit::barsa;
const double distr[3] = { 1.0, 0.0, 0.0 }; // Large enough irrespective of #phases.
well_controls_add_new(BHP, target, distr, ctrl);
well_controls_set_current(ctrl, well_controls_get_num(ctrl) - 1);
}
}
void computeTransportSourceSinglePhase(const UnstructuredGrid& grid,
const std::vector<double>& src,
const std::vector<double>& faceflux,
const double inflow_frac,
const Wells* wells,
const std::vector<double>& well_perfrates,
std::vector<double>& transport_src)
{
using namespace Opm;
int nc = grid.number_of_cells;
transport_src.resize(nc);
// Source term and boundary contributions.
for (int c = 0; c < nc; ++c) {
transport_src[c] = 0.0;
transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c];
for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
int f = grid.cell_faces[hf];
const int* f2c = &grid.face_cells[2*f];
double bdy_influx = 0.0;
if (f2c[0] == c && f2c[1] == -1) {
bdy_influx = -faceflux[f];
} else if (f2c[0] == -1 && f2c[1] == c) {
bdy_influx = faceflux[f];
}
if (bdy_influx != 0.0) {
transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx;
}
}
}
// Well contributions.
if (wells) {
const int nw = wells->number_of_wells;
for (int w = 0; w < nw; ++w) {
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
const int perf_cell = wells->well_cells[perf];
double perf_rate = well_perfrates[perf];
if (perf_rate > 0.0) {
// perf_rate is a total inflow rate, we want a water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
perf_rate *= inflow_frac;
}
}
transport_src[perf_cell] += perf_rate;
}
}
}
}
} // anon namespace } // anon namespace
@ -102,95 +162,34 @@ try
parameter::ParameterGroup param(argc, argv, false); parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl; std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that. // Read the deck.
bool use_deck = param.has("deck_filename"); std::string deck_filename = param.get<std::string>("deck_filename");
Opm::DeckConstPtr deck; Parser parser;
std::unique_ptr<GridManager> grid; DeckConstPtr deck = parser.parseFile(deck_filename);
std::unique_ptr<IncompPropertiesInterface> props; EclipseStateConstPtr eclipseState = std::make_shared<EclipseState>(deck);
std::unique_ptr<Opm::WellsManager> wells;
TwophaseState state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile(deck_filename);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
// Grid init // Grid init
grid.reset(new GridManager(deck)); GridManager grid_manager(deck);
// Rock and fluid init const UnstructuredGrid& grid = *grid_manager.c_grid();
props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid())); // Rock and fluid init
// Wells init. IncompPropertiesSinglePhase props(deck, eclipseState, grid);
wells.reset(new Opm::WellsManager(eclipseState , 0 , *grid->c_grid(), props->permeability())); // Wells init.
// Gravity. WellsManager wells_manager(eclipseState , 0, grid, props.permeability());
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
} 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 GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Wells init.
wells.reset(new Opm::WellsManager());
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
}
// Warn if gravity but no density difference. std::shared_ptr<Wells> my_wells(clone_wells(wells_manager.c_wells()), destroy_wells);
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); setBhpWells(*my_wells);
if (use_gravity) { const Wells& wells = *my_wells;
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
const double *grav = use_gravity ? &gravity[0] : 0;
// Initialising src // Pore volume.
std::vector<double> porevol; std::vector<double> porevol;
computePorevolume(*grid->c_grid(), props->porosity(), porevol); computePorevolume(grid, props.porosity(), porevol);
int num_cells = grid->c_grid()->number_of_cells; int num_cells = grid.number_of_cells;
std::vector<double> src(num_cells, 0.0);
if (use_deck) {
// Do nothing, wells will be the driving force, not source terms.
} else {
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/unit::day;
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
// Boundary conditions.
FlowBCManager bcs;
if (param.getDefault("use_pside", false)) {
int pside = param.get<int>("pside");
double pside_pressure = param.get<double>("pside_pressure");
bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
}
// Linear solver. // Linear solver.
LinearSolverFactory linsolver(param); LinearSolverFactory linsolver(param);
// Pressure solver. // Pressure solver.
Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver, Opm::IncompTpfaSinglePhase psolver(grid, props, linsolver, wells);
0.0, 0.0, 0,
grav, wells->c_wells(), src, bcs.c_bcs());
// Choice of tof solver. // Choice of tof solver.
bool use_dg = param.getDefault("use_dg", false); bool use_dg = param.getDefault("use_dg", false);
@ -198,12 +197,20 @@ try
// Need to initialize dg solver here, since it uses parameters now. // Need to initialize dg solver here, since it uses parameters now.
std::unique_ptr<Opm::TofDiscGalReorder> dg_solver; std::unique_ptr<Opm::TofDiscGalReorder> dg_solver;
if (use_dg) { if (use_dg) {
dg_solver.reset(new Opm::TofDiscGalReorder(*grid->c_grid(), param)); dg_solver.reset(new Opm::TofDiscGalReorder(grid, param));
} else { } else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false); use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
} }
bool compute_tracer = param.getDefault("compute_tracer", false); bool compute_tracer = param.getDefault("compute_tracer", false);
bool start_from_injectors = true;
std::string trace_start = param.getDefault<std::string>("trace_start", "Injectors");
if (trace_start == "Producers") {
start_from_injectors = false;
} else if (trace_start != "Injectors") {
OPM_THROW(std::runtime_error, "Unknown trace_start specification (allowed is Injectors or Producers): " << trace_start);
}
// Write parameters used for later reference. // Write parameters used for later reference.
bool output = param.getDefault("output", true); bool output = param.getDefault("output", true);
std::ofstream epoch_os; std::ofstream epoch_os;
@ -218,19 +225,9 @@ try
catch (...) { catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
} }
std::string filename = output_dir + "/epoch_timing.param";
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
// open file to clean it. The file is appended to in SimulatorTwophase
filename = output_dir + "/step_timing.param";
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
step_os.close();
param.writeParam(output_dir + "/simulation.param"); param.writeParam(output_dir + "/simulation.param");
} }
// Init wells.
Opm::WellState well_state;
well_state.init(wells->c_wells(), state);
// Check if we have misspelled anything // Check if we have misspelled anything
warnIfUnusedParams(param); warnIfUnusedParams(param);
@ -244,17 +241,29 @@ try
std::cout << "\n\n================ Starting main solvers ===============" << std::endl; std::cout << "\n\n================ Starting main solvers ===============" << std::endl;
// Solve pressure. // Solve pressure.
std::vector<double> press;
std::vector<double> flux;
std::vector<double> bhp;
std::vector<double> wellrates;
pressure_timer.start(); pressure_timer.start();
psolver.solve(1.0, state, well_state); psolver.solve(press, flux, bhp, wellrates);
pressure_timer.stop(); pressure_timer.stop();
double pt = pressure_timer.secsSinceStart(); double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt; ptime += pt;
// Turn direction of flux if starting from producers.
if (!start_from_injectors) {
for (auto it = flux.begin(); it != flux.end(); ++it) {
(*it) = -(*it);
}
}
// Process transport sources (to include bdy terms and well flows). // Process transport sources (to include bdy terms and well flows).
std::vector<double> src(num_cells, 0.0);
std::vector<double> transport_src; std::vector<double> transport_src;
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, computeTransportSourceSinglePhase(grid, src, flux, 1.0,
wells->c_wells(), well_state.perfRates(), transport_src); &wells, wellrates, transport_src);
// Solve time-of-flight. // Solve time-of-flight.
transport_timer.start(); transport_timer.start();
@ -262,20 +271,20 @@ try
std::vector<double> tracer; std::vector<double> tracer;
Opm::SparseTable<int> tracerheads; Opm::SparseTable<int> tracerheads;
if (compute_tracer) { if (compute_tracer) {
buildTracerheadsFromWells(wells->c_wells(), tracerheads); buildTracerheadsFromWells(wells, start_from_injectors, tracerheads);
} }
if (use_dg) { if (use_dg) {
if (compute_tracer) { if (compute_tracer) {
dg_solver->solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tracerheads, tof, tracer); dg_solver->solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else { } else {
dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); dg_solver->solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
} }
} else { } else {
Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); Opm::TofReorder tofsolver(grid, use_multidim_upwind);
if (compute_tracer) { if (compute_tracer) {
tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tracerheads, tof, tracer); tofsolver.solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else { } else {
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); tofsolver.solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
} }
} }
transport_timer.stop(); transport_timer.stop();

155
opm/core/pressure/IncompTpfaSinglePhase.cpp Normal file
View File

@ -0,0 +1,155 @@
/*
Copyright 2015 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 "config.h"
#include <opm/core/pressure/IncompTpfaSinglePhase.hpp>
#include <opm/core/props/IncompPropertiesSinglePhase.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
// #include <opm/core/pressure/mimetic/mimetic.h>
// #include <opm/core/pressure/flow_bc.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/sparse_sys.h>
// #include <opm/core/simulator/TwophaseState.hpp>
// #include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
// #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/wells.h>
// #include <iostream>
// #include <iomanip>
// #include <cmath>
// #include <algorithm>
namespace Opm
{
/// Construct solver for incompressible case.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver to use.
/// \param[in] wells The wells used as driving forces.
IncompTpfaSinglePhase::IncompTpfaSinglePhase(const UnstructuredGrid& grid,
const IncompPropertiesSinglePhase& props,
const LinearSolverInterface& linsolver,
const Wells& wells)
: grid_(grid),
props_(props),
linsolver_(linsolver),
wells_(wells),
htrans_(grid.cell_facepos[ grid.number_of_cells ]),
trans_ (grid.number_of_faces),
zeros_(grid.cell_facepos[ grid.number_of_cells ])
{
computeStaticData();
}
/// Destructor.
IncompTpfaSinglePhase::~IncompTpfaSinglePhase()
{
ifs_tpfa_destroy(h_);
}
/// Solve the pressure equation.
void IncompTpfaSinglePhase::solve(std::vector<double>& press,
std::vector<double>& flux,
std::vector<double>& bhp,
std::vector<double>& wellrates)
{
// Set up properties.
computePerSolveDynamicData();
// Assemble.
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
int ok = ifs_tpfa_assemble(gg, &forces_, trans_.data(), zeros_.data(), h_);
if (!ok) {
OPM_THROW(std::runtime_error, "Failed assembling pressure system.");
}
// Solve.
linsolver_.solve(h_->A, h_->b, h_->x);
// Obtain solution.
press.resize(grid_.number_of_cells);
flux.resize(grid_.number_of_faces);
wellrates.resize(wells_.well_connpos[ wells_.number_of_wells ]);
bhp.resize(wells_.number_of_wells);
ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
soln.cell_press = press.data();
soln.face_flux = flux.data();
soln.well_press = bhp.data();
soln.well_flux = wellrates.data();
ifs_tpfa_press_flux(gg, &forces_, &trans_[0], h_, &soln);
}
/// Compute data that never changes (after construction).
void IncompTpfaSinglePhase::computeStaticData()
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_htrans_compute(gg, props_.permeability(), &htrans_[0]);
h_ = ifs_tpfa_construct(gg, const_cast<struct Wells*>(&wells_));
}
/// Compute per-solve dynamic properties.
void IncompTpfaSinglePhase::computePerSolveDynamicData()
{
// Computed here:
//
// std::vector<double> totmob_;
// std::vector<double> trans_;
// ifs_tpfa_forces forces_;
// totmob_
totmob_.clear();
totmob_.resize(grid_.number_of_cells, 1.0/(*props_.viscosity()));
// trans_
tpfa_eff_trans_compute(const_cast<UnstructuredGrid*>(&grid_), totmob_.data(), htrans_.data(), trans_.data());
// forces_
forces_.src = NULL;
forces_.bc = NULL;
forces_.W = &wells_;
forces_.totmob = totmob_.data();
forces_.wdp = zeros_.data();
}
} // namespace Opm

88
opm/core/pressure/IncompTpfaSinglePhase.hpp Normal file
View File

@ -0,0 +1,88 @@
/*
Copyright 2015 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_INCOMPTPFASINGLEPHASE_HEADER_INCLUDED
#define OPM_INCOMPTPFASINGLEPHASE_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <vector>
struct UnstructuredGrid;
struct Wells;
namespace Opm
{
class IncompPropertiesSinglePhase;
class LinearSolverInterface;
/// Encapsulating a tpfa pressure solver for the incompressible-fluid case.
/// Supports gravity, wells controlled by bhp or reservoir rates,
/// boundary conditions and simple sources as driving forces.
/// Rock compressibility can be included, and necessary nonlinear
/// iterations are handled.
/// Below we use the shortcuts D for the number of dimensions, N
/// for the number of cells and F for the number of faces.
class IncompTpfaSinglePhase
{
public:
/// Construct solver for incompressible case.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver to use.
/// \param[in] wells The wells used as driving forces.
IncompTpfaSinglePhase(const UnstructuredGrid& grid,
const IncompPropertiesSinglePhase& props,
const LinearSolverInterface& linsolver,
const Wells& wells);
/// Destructor.
~IncompTpfaSinglePhase();
/// Solve the pressure equation.
void solve(std::vector<double>& press,
std::vector<double>& flux,
std::vector<double>& bhp,
std::vector<double>& wellrates);
private:
// Helper functions.
void computeStaticData();
void computePerSolveDynamicData();
protected:
// ------ Data that will remain unmodified after construction. ------
const UnstructuredGrid& grid_;
const IncompPropertiesSinglePhase& props_;
const LinearSolverInterface& linsolver_;
const Wells& wells_;
std::vector<double> htrans_;
std::vector<double> trans_ ;
std::vector<double> zeros_;
std::vector<double> totmob_;
struct ifs_tpfa_forces forces_;
// ------ Internal data for the ifs_tpfa solver. ------
struct ifs_tpfa_data* h_;
};
} // namespace Opm
#endif // OPM_INCOMPTPFASINGLEPHASE_HEADER_INCLUDED

178
opm/core/props/IncompPropertiesSinglePhase.cpp Normal file
View File

@ -0,0 +1,178 @@
/*
Copyright 2015 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 "config.h"
#include <opm/core/props/IncompPropertiesSinglePhase.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
// #include <iostream>
namespace Opm
{
IncompPropertiesSinglePhase::IncompPropertiesSinglePhase(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
const UnstructuredGrid& grid)
{
rock_.init(eclState, grid.number_of_cells, grid.global_cell, grid.cartdims);
if (deck->hasKeyword("DENSITY")) {
Opm::DeckRecordConstPtr densityRecord = deck->getKeyword("DENSITY")->getRecord(0);
surface_density_ = densityRecord->getItem("OIL")->getSIDouble(0);
} else {
surface_density_ = 1000.0;
OPM_MESSAGE("Input is missing DENSITY -- using a standard density of "
<< surface_density_ << ".\n");
}
// This will be modified if we have a PVCDO specification.
reservoir_density_ = surface_density_;
if (deck->hasKeyword("PVCDO")) {
Opm::DeckRecordConstPtr pvcdoRecord = deck->getKeyword("PVCDO")->getRecord(0);
if (pvcdoRecord->getItem("OIL_COMPRESSIBILITY")->getSIDouble(0) != 0.0 ||
pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0) != 0.0) {
OPM_MESSAGE("Compressibility effects in PVCDO are ignored.");
}
reservoir_density_ /= pvcdoRecord->getItem("OIL_VOL_FACTOR")->getSIDouble(0);
viscosity_ = pvcdoRecord->getItem("OIL_VISCOSITY")->getSIDouble(0);
} else {
viscosity_ = 1.0 * prefix::centi*unit::Poise;
OPM_MESSAGE("Input is missing PVCDO -- using a standard viscosity of "
<< viscosity_ << " and reservoir density equal to surface density.\n");
}
}
IncompPropertiesSinglePhase::~IncompPropertiesSinglePhase()
{
}
/// \return D, the number of spatial dimensions.
int IncompPropertiesSinglePhase::numDimensions() const
{
return rock_.numDimensions();
}
/// \return N, the number of cells.
int IncompPropertiesSinglePhase::numCells() const
{
return rock_.numCells();
}
/// \return Array of N porosity values.
const double* IncompPropertiesSinglePhase::porosity() const
{
return rock_.porosity();
}
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* IncompPropertiesSinglePhase::permeability() const
{
return rock_.permeability();
}
// ---- Fluid interface ----
/// \return P, the number of phases (also the number of components).
int IncompPropertiesSinglePhase::numPhases() const
{
return 1;
}
/// \return Array of P viscosity values.
const double* IncompPropertiesSinglePhase::viscosity() const
{
return &viscosity_;
}
/// \return Array of P density values.
const double* IncompPropertiesSinglePhase::density() const
{
return &reservoir_density_;
}
/// \return Array of P density values.
const double* IncompPropertiesSinglePhase::surfaceDensity() const
{
return &surface_density_;
}
/// Relative permeability. Always returns 1 (and 0 for derivatives).
/// \param[in] n Number of data points.
/// \param[in] s Array of n saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of n relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of n relperm derivative values,
/// array must be valid before calling.
void IncompPropertiesSinglePhase::relperm(const int n,
const double* /* s */,
const int* /* cells */,
double* kr,
double* dkrds) const
{
std::fill(kr, kr + n, 1.0);
if (dkrds) {
std::fill(dkrds, dkrds + n, 0.0);
}
}
/// Capillary pressure. Always returns zero.
/// \param[in] n Number of data points.
/// \param[in] s Array of n saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of n capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of n derivative values,
/// array must be valid before calling.
void IncompPropertiesSinglePhase::capPress(const int n,
const double* /* s */,
const int* /* cells */,
double* pc,
double* dpcds) const
{
std::fill(pc, pc + n, 0.0);
if (dpcds) {
std::fill(dpcds, dpcds + n, 0.0);
}
}
/// Obtain the range of allowable saturation values.
/// Saturation range is just the point 1 for this class
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of n minimum s values, array must be valid before calling.
/// \param[out] smax Array of n maximum s values, array must be valid before calling.
void IncompPropertiesSinglePhase::satRange(const int n,
const int* /* cells */,
double* smin,
double* smax) const
{
std::fill(smin, smin + n, 1.0);
std::fill(smax, smax + n, 1.0);
}
} // namespace Opm

146
opm/core/props/IncompPropertiesSinglePhase.hpp Normal file
View File

@ -0,0 +1,146 @@
/*
Copyright 2015 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_INCOMPPROPERTIESSINGLEPHASE_HEADER_INCLUDED
#define OPM_INCOMPPROPERTIESSINGLEPHASE_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockFromDeck.hpp>
struct UnstructuredGrid;
namespace Opm
{
/// Concrete class implementing the incompressible property
/// interface for a simplified single-phase setting, reading all
/// data and properties from eclipse deck input. The oil phase
/// properties are used where applicable and available.
///
/// Supports variable number of spatial dimensions, called D.
/// Supports a single phase only.
/// In general, when arguments call for n values of some vector or
/// matrix property, such as saturation, they shall always be
/// ordered cellwise:
/// [s^1_0 s^2_0 s^3_0 s^1_1 s^2_2 ... ]
/// in which s^i_j denotes saturation of phase i in cell j.
class IncompPropertiesSinglePhase : public IncompPropertiesInterface
{
public:
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param eclState The EclipseState (processed deck) produced by the opm-parser code
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
IncompPropertiesSinglePhase(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
const UnstructuredGrid& grid);
/// Destructor.
virtual ~IncompPropertiesSinglePhase();
// ---- Rock interface ----
/// \return D, the number of spatial dimensions.
virtual int numDimensions() const;
/// \return N, the number of cells.
virtual int numCells() const;
/// \return Array of N porosity values.
virtual const double* porosity() const;
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
virtual const double* permeability() const;
// ---- Fluid interface ----
/// \return P, the number of phases (= 1).
virtual int numPhases() const;
/// \return Array of P (= 1) viscosity values.
virtual const double* viscosity() const;
/// Densities of fluid at reservoir conditions.
/// \return Array of P (= 1) density values.
virtual const double* density() const;
/// Densities of fluid phases at surface conditions.
/// \return Array of P (= 1) density values.
virtual const double* surfaceDensity() const;
/// Relative permeability. Always returns 1 (and 0 for derivatives).
/// \param[in] n Number of data points.
/// \param[in] s Array of n saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] kr Array of n relperm values, array must be valid before calling.
/// \param[out] dkrds If non-null: array of n relperm derivative values,
/// array must be valid before calling.
virtual void relperm(const int n,
const double* s,
const int* cells,
double* kr,
double* dkrds) const;
/// Capillary pressure. Always returns zero.
/// \param[in] n Number of data points.
/// \param[in] s Array of n saturation values.
/// \param[in] cells Array of n cell indices to be associated with the s values.
/// \param[out] pc Array of n capillary pressure values, array must be valid before calling.
/// \param[out] dpcds If non-null: array of n derivative values,
/// array must be valid before calling.
virtual void capPress(const int n,
const double* s,
const int* cells,
double* pc,
double* dpcds) const;
/// Obtain the range of allowable saturation values.
/// Saturation range is just the point 1 for this class
/// \param[in] n Number of data points.
/// \param[in] cells Array of n cell indices.
/// \param[out] smin Array of n minimum s values, array must be valid before calling.
/// \param[out] smax Array of n maximum s values, array must be valid before calling.
virtual void satRange(const int n,
const int* cells,
double* smin,
double* smax) const;
private:
RockFromDeck rock_;
double surface_density_;
double reservoir_density_;
double viscosity_;
};
} // namespace Opm
#endif // OPM_INCOMPPROPERTIESSINGLEPHASE_HEADER_INCLUDED

166
opm/core/tof/FlowDiagnostics.cpp Normal file
View File

@ -0,0 +1,166 @@
/*
Copyright 2015 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/tof/FlowDiagnostics.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <numeric>
namespace Opm
{
/// \brief Compute flow-capacity/storage-capacity based on time-of-flight.
///
/// The F-Phi curve is an analogue to the fractional flow curve in a 1D
/// displacement. It can be used to compute other interesting diagnostic
/// quantities such as the Lorenz coefficient. For a technical description
/// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625).
///
/// \param[in] pv pore volumes of each cell
/// \param[in] ftof forward (time from injector) time-of-flight values for each cell
/// \param[in] rtof reverse (time to producer) time-of-flight values for each cell
/// \return a pair of vectors, the first containing F (flow capacity) the second
/// containing Phi (storage capacity).
std::pair<std::vector<double>, std::vector<double>> computeFandPhi(const std::vector<double>& pv,
const std::vector<double>& ftof,
const std::vector<double>& rtof)
{
if (pv.size() != ftof.size() || pv.size() != rtof.size()) {
OPM_THROW(std::runtime_error, "computeFandPhi(): Input vectors must have same size.");
}
// Sort according to total travel time.
const int n = pv.size();
typedef std::pair<double, double> D2;
std::vector<D2> time_and_pv(n);
for (int ii = 0; ii < n; ++ii) {
time_and_pv[ii].first = ftof[ii] + rtof[ii]; // Total travel time.
time_and_pv[ii].second = pv[ii];
}
std::sort(time_and_pv.begin(), time_and_pv.end());
// Compute Phi.
std::vector<double> Phi(n + 1);
Phi[0] = 0.0;
for (int ii = 0; ii < n; ++ii) {
Phi[ii+1] = time_and_pv[ii].second;
}
std::partial_sum(Phi.begin(), Phi.end(), Phi.begin());
const double vt = Phi.back(); // Total pore volume.
for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop.
Phi[ii] /= vt; // Normalize Phi.
}
// Compute F.
std::vector<double> F(n + 1);
F[0] = 0.0;
for (int ii = 0; ii < n; ++ii) {
F[ii+1] = time_and_pv[ii].second / time_and_pv[ii].first;
}
std::partial_sum(F.begin(), F.end(), F.begin());
const double ft = F.back(); // Total flux.
for (int ii = 1; ii < n+1; ++ii) { // Note limits of loop.
F[ii] /= ft; // Normalize Phi.
}
return std::make_pair(F, Phi);
}
/// \brief Compute the Lorenz coefficient based on the F-Phi curve.
///
/// The Lorenz coefficient is a measure of heterogeneity. It is equal
/// to twice the area between the F-Phi curve and the F = Phi line.
/// The coefficient can vary from zero to one. If the coefficient is
/// zero (so the F-Phi curve is a straight line) we have perfect
/// piston-like displacement while a coefficient of one indicates
/// infinitely heterogenous displacement (essentially no sweep).
///
/// Note: The coefficient is analogous to the Gini coefficient of
/// economic theory, where the name Lorenz curve is applied to
/// what we call the F-Phi curve.
///
/// \param[in] flowcap flow capacity (F) as from computeFandPhi()
/// \param[in] storagecap storage capacity (Phi) as from computeFandPhi()
/// \return the Lorenz coefficient
double computeLorenz(const std::vector<double>& flowcap,
const std::vector<double>& storagecap)
{
if (flowcap.size() != storagecap.size()) {
OPM_THROW(std::runtime_error, "computeLorenz(): Input vectors must have same size.");
}
double integral = 0.0;
// Trapezoid quadrature of the curve F(Phi).
const int num_intervals = flowcap.size() - 1;
for (int ii = 0; ii < num_intervals; ++ii) {
const double len = storagecap[ii+1] - storagecap[ii];
integral += (flowcap[ii] + flowcap[ii+1]) * len / 2.0;
}
return 2.0 * (integral - 0.5);
}
/// \brief Compute sweep efficiency versus dimensionless time (PVI).
///
/// The sweep efficiency is analogue to 1D displacement using the
/// F-Phi curve as flux function.
///
/// \param[in] flowcap flow capacity (F) as from computeFandPhi()
/// \param[in] storagecap storage capacity (Phi) as from computeFandPhi()
/// \return a pair of vectors, the first containing Ev (sweep efficiency)
/// the second containing tD (dimensionless time).
std::pair<std::vector<double>, std::vector<double>> computeSweep(const std::vector<double>& flowcap,
const std::vector<double>& storagecap)
{
if (flowcap.size() != storagecap.size()) {
OPM_THROW(std::runtime_error, "computeSweep(): Input vectors must have same size.");
}
// Compute tD and Ev simultaneously,
// skipping identical Phi data points.
const int n = flowcap.size();
std::vector<double> Ev;
std::vector<double> tD;
tD.reserve(n);
Ev.reserve(n);
tD.push_back(0.0);
Ev.push_back(0.0);
for (int ii = 1; ii < n; ++ii) { // Note loop limits.
const double fd = flowcap[ii] - flowcap[ii-1];
const double sd = storagecap[ii] - storagecap[ii-1];
if (fd != 0.0) {
tD.push_back(sd/fd);
Ev.push_back(storagecap[ii] + (1.0 - flowcap[ii]) * tD.back());
}
}
return std::make_pair(Ev, tD);
}
} // namespace Opm

84
opm/core/tof/FlowDiagnostics.hpp Normal file
View File

@ -0,0 +1,84 @@
/*
Copyright 2015 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_FLOWDIAGNOSTICS_HEADER_INCLUDED
#define OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED
#include <vector>
#include <utility>
namespace Opm
{
/// \brief Compute flow-capacity/storage-capacity based on time-of-flight.
///
/// The F-Phi curve is an analogue to the fractional flow curve in a 1D
/// displacement. It can be used to compute other interesting diagnostic
/// quantities such as the Lorenz coefficient. For a technical description
/// see Shavali et al. (SPE 146446), Shook and Mitchell (SPE 124625).
///
/// \param[in] pv pore volumes of each cell
/// \param[in] ftof forward (time from injector) time-of-flight values for each cell
/// \param[in] rtof reverse (time to producer) time-of-flight values for each cell
/// \return a pair of vectors, the first containing F (flow capacity) the second
/// containing Phi (storage capacity).
std::pair<std::vector<double>, std::vector<double>>
computeFandPhi(const std::vector<double>& pv,
const std::vector<double>& ftof,
const std::vector<double>& rtof);
/// \brief Compute the Lorenz coefficient based on the F-Phi curve.
///
/// The Lorenz coefficient is a measure of heterogeneity. It is equal
/// to twice the area between the F-Phi curve and the F = Phi line.
/// The coefficient can vary from zero to one. If the coefficient is
/// zero (so the F-Phi curve is a straight line) we have perfect
/// piston-like displacement while a coefficient of one indicates
/// infinitely heterogenous displacement (essentially no sweep).
///
/// Note: The coefficient is analogous to the Gini coefficient of
/// economic theory, where the name Lorenz curve is applied to
/// what we call the F-Phi curve.
///
/// \param[in] flowcap flow capacity (F) as from computeFandPhi()
/// \param[in] storagecap storage capacity (Phi) as from computeFandPhi()
/// \return the Lorenz coefficient
double computeLorenz(const std::vector<double>& flowcap,
const std::vector<double>& storagecap);
/// \brief Compute sweep efficiency versus dimensionless time (PVI).
///
/// The sweep efficiency is analogue to 1D displacement using the
/// F-Phi curve as flux function.
///
/// \param[in] flowcap flow capacity (F) as from computeFandPhi()
/// \param[in] storagecap storage capacity (Phi) as from computeFandPhi()
/// \return a pair of vectors, the first containing Ev (sweep efficiency)
/// the second containing tD (dimensionless time).
std::pair<std::vector<double>, std::vector<double>>
computeSweep(const std::vector<double>& flowcap,
const std::vector<double>& storagecap);
} // namespace Opm
#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED

198
tests/test_flowdiagnostics.cpp Normal file
View File

@ -0,0 +1,198 @@
/*
Copyright 2015 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 <config.h>
#if defined(HAVE_DYNAMIC_BOOST_TEST)
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE FlowDiagnosticsTests
#include <boost/test/unit_test.hpp>
#include <opm/core/tof/FlowDiagnostics.hpp>
const std::vector<double> pv(16, 18750.0);
const std::vector<double> ftof = {
5.399999999999992e+04,
1.139999999999997e+05,
2.819999999999993e+05,
8.220000000000012e+05,
1.139999999999998e+05,
1.774285714285711e+05,
3.160150375939849e+05,
8.156820645778908e+05,
2.819999999999994e+05,
3.160150375939841e+05,
3.935500938781204e+05,
7.765612369042073e+05,
8.220000000000000e+05,
8.156820645778894e+05,
7.765612369042063e+05,
8.218220225906991e+05
};
const std::vector<double> rtof = {
8.218220225906990e+05,
7.765612369042046e+05,
8.156820645778881e+05,
8.219999999999976e+05,
7.765612369042051e+05,
3.935500938781204e+05,
3.160150375939846e+05,
2.820000000000001e+05,
8.156820645778885e+05,
3.160150375939850e+05,
1.774285714285714e+05,
1.140000000000000e+05,
8.219999999999980e+05,
2.819999999999998e+05,
1.140000000000000e+05,
5.400000000000003e+04
};
const std::vector<double> F = {
0,
9.568875799840706e-02,
1.913775159968141e-01,
2.778231480508526e-01,
3.642687801048911e-01,
4.266515906731506e-01,
4.890344012414101e-01,
5.503847464649610e-01,
6.117350916885119e-01,
6.730854369120627e-01,
7.344357821356134e-01,
7.842099754925904e-01,
8.339841688495674e-01,
8.837583622065442e-01,
9.335325555635212e-01,
9.667662777817606e-01,
1.000000000000000e+00
};
const std::vector<double> Phi = {
0,
6.250000000000000e-02,
1.250000000000000e-01,
1.875000000000000e-01,
2.500000000000000e-01,
3.125000000000000e-01,
3.750000000000000e-01,
4.375000000000000e-01,
5.000000000000000e-01,
5.625000000000000e-01,
6.250000000000000e-01,
6.875000000000000e-01,
7.500000000000000e-01,
8.125000000000000e-01,
8.750000000000000e-01,
9.375000000000000e-01,
1.000000000000000e+00
};
const std::vector<double> Ev = {
0,
6.531592770912591e-01,
6.531592770912593e-01,
7.096322601771997e-01,
7.096322601772002e-01,
8.869254748464411e-01,
8.869254748464422e-01,
8.955406718746977e-01,
8.955406718746983e-01,
8.955406718746991e-01,
8.955406718746991e-01,
9.584612275378565e-01,
9.584612275378565e-01,
9.584612275378569e-01,
9.584612275378566e-01,
1.000000000000000e+00,
1.000000000000000e+00
};
const std::vector<double> tD = {
0,
6.531592770912591e-01,
6.531592770912593e-01,
7.229977792392133e-01,
7.229977792392139e-01,
1.001878553253259e+00,
1.001878553253261e+00,
1.018739173712224e+00,
1.018739173712226e+00,
1.018739173712227e+00,
1.018739173712227e+00,
1.255670776053656e+00,
1.255670776053656e+00,
1.255670776053659e+00,
1.255670776053656e+00,
1.880619919417231e+00,
1.880619919417231e+00
};
std::vector<double> wrong_length(7, 0.0);
using namespace Opm;
template <class C>
void compareCollections(const C& c1, const C& c2, const double tolerance = 1e-11)
{
BOOST_REQUIRE(c1.size() == c2.size());
auto c1it = c1.begin();
auto c2it = c2.begin();
for (; c1it != c1.end(); ++c1it, ++c2it) {
BOOST_CHECK_CLOSE(*c1it, *c2it, tolerance);
}
}
BOOST_AUTO_TEST_CASE(FandPhi)
{
BOOST_CHECK_THROW(computeFandPhi(pv, ftof, wrong_length), std::runtime_error);
auto FPhi = computeFandPhi(pv, ftof, rtof);
compareCollections(FPhi.first, F);
compareCollections(FPhi.second, Phi);
}
BOOST_AUTO_TEST_CASE(Lorenz)
{
BOOST_CHECK_THROW(computeLorenz(F, wrong_length), std::runtime_error);
const double Lc = computeLorenz(F, Phi);
BOOST_CHECK_CLOSE(Lc, 1.645920738950826e-01, 1e-11);
}
BOOST_AUTO_TEST_CASE(Sweep)
{
BOOST_CHECK_THROW(computeSweep(F, wrong_length), std::runtime_error);
auto et = computeSweep(F, Phi);
compareCollections(et.first, Ev);
compareCollections(et.second, tD);
}