diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 149f73f9..debb3de6 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -55,6 +55,7 @@ list (APPEND MAIN_SOURCE_FILES opm/core/pressure/CompressibleTpfa.cpp opm/core/pressure/FlowBCManager.cpp opm/core/pressure/IncompTpfa.cpp + opm/core/pressure/IncompTpfaSinglePhase.cpp opm/core/pressure/cfsh.c opm/core/pressure/flow_bc.c opm/core/pressure/fsh.c @@ -83,6 +84,7 @@ list (APPEND MAIN_SOURCE_FILES opm/core/props/BlackoilPropertiesFromDeck.cpp opm/core/props/IncompPropertiesBasic.cpp opm/core/props/IncompPropertiesFromDeck.cpp + opm/core/props/IncompPropertiesSinglePhase.cpp opm/core/props/pvt/BlackoilPvtProperties.cpp opm/core/props/pvt/PvtPropertiesBasic.cpp opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp @@ -110,6 +112,7 @@ list (APPEND MAIN_SOURCE_FILES opm/core/simulator/SimulatorTimer.cpp opm/core/tof/AnisotropicEikonal.cpp opm/core/tof/DGBasis.cpp + opm/core/tof/FlowDiagnostics.cpp opm/core/tof/TofReorder.cpp opm/core/tof/TofDiscGalReorder.cpp opm/core/transport/TransportSolverTwophaseInterface.cpp @@ -162,6 +165,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_ug.cpp tests/test_cubic.cpp tests/test_event.cpp + tests/test_flowdiagnostics.cpp tests/test_nonuniformtablelinear.cpp tests/test_sparsevector.cpp tests/test_sparsetable.cpp @@ -302,6 +306,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/pressure/CompressibleTpfa.hpp opm/core/pressure/FlowBCManager.hpp opm/core/pressure/IncompTpfa.hpp + opm/core/pressure/IncompTpfaSinglePhase.hpp opm/core/pressure/flow_bc.h opm/core/pressure/fsh.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/IncompPropertiesShadow.hpp opm/core/props/IncompPropertiesShadow_impl.hpp + opm/core/props/IncompPropertiesSinglePhase.hpp opm/core/props/phaseUsageFromDeck.hpp opm/core/props/pvt/BlackoilPvtProperties.hpp opm/core/props/pvt/PvtPropertiesBasic.hpp @@ -379,6 +385,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/simulator/initStateEquil_impl.hpp opm/core/tof/AnisotropicEikonal.hpp opm/core/tof/DGBasis.hpp + opm/core/tof/FlowDiagnostics.hpp opm/core/tof/TofReorder.hpp opm/core/tof/TofDiscGalReorder.hpp opm/core/transport/TransportSolverTwophaseInterface.hpp diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 087e10b9..7dd20f8b 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -23,8 +23,6 @@ #include "config.h" #endif // HAVE_CONFIG_H -#include - #include #include #include @@ -32,18 +30,15 @@ #include #include #include +#include #include #include -#include -#include +#include #include -#include -#include -#include -#include +#include #include #include @@ -70,23 +65,88 @@ namespace } } - void buildTracerheadsFromWells(const Wells* wells, + void buildTracerheadsFromWells(const Wells& wells, + const bool trace_injectors, Opm::SparseTable& tracerheads) { - if (wells == 0) { - return; - } 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) { - if (wells->type[w] != INJECTOR) { + if (wells.type[w] != wanted_type) { continue; } - tracerheads.appendRow(wells->well_cells + wells->well_connpos[w], - wells->well_cells + wells->well_connpos[w + 1]); + tracerheads.appendRow(wells.well_cells + wells.well_connpos[w], + 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& src, + const std::vector& faceflux, + const double inflow_frac, + const Wells* wells, + const std::vector& well_perfrates, + std::vector& 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 @@ -102,95 +162,34 @@ try parameter::ParameterGroup param(argc, argv, false); std::cout << "--------------- Reading parameters ---------------" << std::endl; - // If we have a "deck_filename", grid and props will be read from that. - bool use_deck = param.has("deck_filename"); - Opm::DeckConstPtr deck; - std::unique_ptr grid; - std::unique_ptr props; - std::unique_ptr 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("deck_filename"); - Opm::ParserPtr parser(new Opm::Parser()); - deck = parser->parseFile(deck_filename); - Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); + // Read the deck. + std::string deck_filename = param.get("deck_filename"); + Parser parser; + DeckConstPtr deck = parser.parseFile(deck_filename); + EclipseStateConstPtr eclipseState = std::make_shared(deck); - // Grid init - grid.reset(new GridManager(deck)); - // Rock and fluid init - props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid())); - // Wells init. - wells.reset(new Opm::WellsManager(eclipseState , 0 , *grid->c_grid(), props->permeability())); - // Gravity. - 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); - } + // Grid init + GridManager grid_manager(deck); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + // Rock and fluid init + IncompPropertiesSinglePhase props(deck, eclipseState, grid); + // Wells init. + WellsManager wells_manager(eclipseState , 0, grid, props.permeability()); - // Warn if gravity but no density difference. - bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); - if (use_gravity) { - 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; + std::shared_ptr my_wells(clone_wells(wells_manager.c_wells()), destroy_wells); + setBhpWells(*my_wells); + const Wells& wells = *my_wells; - // Initialising src + // Pore volume. std::vector porevol; - computePorevolume(*grid->c_grid(), props->porosity(), porevol); - int num_cells = grid->c_grid()->number_of_cells; - std::vector 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("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("pside"); - double pside_pressure = param.get("pside_pressure"); - bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure); - } + computePorevolume(grid, props.porosity(), porevol); + int num_cells = grid.number_of_cells; // Linear solver. LinearSolverFactory linsolver(param); // Pressure solver. - Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver, - 0.0, 0.0, 0, - grav, wells->c_wells(), src, bcs.c_bcs()); + Opm::IncompTpfaSinglePhase psolver(grid, props, linsolver, wells); // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); @@ -198,12 +197,20 @@ try // Need to initialize dg solver here, since it uses parameters now. std::unique_ptr dg_solver; if (use_dg) { - dg_solver.reset(new Opm::TofDiscGalReorder(*grid->c_grid(), param)); + dg_solver.reset(new Opm::TofDiscGalReorder(grid, param)); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } bool compute_tracer = param.getDefault("compute_tracer", false); + bool start_from_injectors = true; + std::string trace_start = param.getDefault("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. bool output = param.getDefault("output", true); std::ofstream epoch_os; @@ -218,19 +225,9 @@ try catch (...) { 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"); } - // Init wells. - Opm::WellState well_state; - well_state.init(wells->c_wells(), state); - // Check if we have misspelled anything warnIfUnusedParams(param); @@ -244,17 +241,29 @@ try std::cout << "\n\n================ Starting main solvers ===============" << std::endl; // Solve pressure. + std::vector press; + std::vector flux; + std::vector bhp; + std::vector wellrates; pressure_timer.start(); - psolver.solve(1.0, state, well_state); + psolver.solve(press, flux, bhp, wellrates); pressure_timer.stop(); double pt = pressure_timer.secsSinceStart(); std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; 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). + std::vector src(num_cells, 0.0); std::vector transport_src; - Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, - wells->c_wells(), well_state.perfRates(), transport_src); + computeTransportSourceSinglePhase(grid, src, flux, 1.0, + &wells, wellrates, transport_src); // Solve time-of-flight. transport_timer.start(); @@ -262,20 +271,20 @@ try std::vector tracer; Opm::SparseTable tracerheads; if (compute_tracer) { - buildTracerheadsFromWells(wells->c_wells(), tracerheads); + buildTracerheadsFromWells(wells, start_from_injectors, tracerheads); } if (use_dg) { 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 { - dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + dg_solver->solveTof(flux.data(), porevol.data(), transport_src.data(), tof); } } else { - Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); + Opm::TofReorder tofsolver(grid, use_multidim_upwind); 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 { - tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + tofsolver.solveTof(flux.data(), porevol.data(), transport_src.data(), tof); } } transport_timer.stop(); diff --git a/opm/core/pressure/IncompTpfaSinglePhase.cpp b/opm/core/pressure/IncompTpfaSinglePhase.cpp new file mode 100644 index 00000000..b39778c6 --- /dev/null +++ b/opm/core/pressure/IncompTpfaSinglePhase.cpp @@ -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 . +*/ + + +#include "config.h" +#include + +#include +#include +#include +// #include +// #include +#include +#include +// #include +// #include +#include +// #include +#include +// #include +// #include +// #include +// #include + +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& press, + std::vector& flux, + std::vector& bhp, + std::vector& wellrates) + { + // Set up properties. + computePerSolveDynamicData(); + + // Assemble. + UnstructuredGrid* gg = const_cast(&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(&grid_); + tpfa_htrans_compute(gg, props_.permeability(), &htrans_[0]); + h_ = ifs_tpfa_construct(gg, const_cast(&wells_)); + } + + + + + + + /// Compute per-solve dynamic properties. + void IncompTpfaSinglePhase::computePerSolveDynamicData() + { + // Computed here: + // + // std::vector totmob_; + // std::vector 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(&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 diff --git a/opm/core/pressure/IncompTpfaSinglePhase.hpp b/opm/core/pressure/IncompTpfaSinglePhase.hpp new file mode 100644 index 00000000..617a093c --- /dev/null +++ b/opm/core/pressure/IncompTpfaSinglePhase.hpp @@ -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 . +*/ + +#ifndef OPM_INCOMPTPFASINGLEPHASE_HEADER_INCLUDED +#define OPM_INCOMPTPFASINGLEPHASE_HEADER_INCLUDED + + +#include +#include + +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& press, + std::vector& flux, + std::vector& bhp, + std::vector& 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 htrans_; + std::vector trans_ ; + std::vector zeros_; + std::vector 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 diff --git a/opm/core/props/IncompPropertiesSinglePhase.cpp b/opm/core/props/IncompPropertiesSinglePhase.cpp new file mode 100644 index 00000000..e56ed32e --- /dev/null +++ b/opm/core/props/IncompPropertiesSinglePhase.cpp @@ -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 . +*/ + + +#include "config.h" + +#include +#include +#include +#include +// #include + +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 + diff --git a/opm/core/props/IncompPropertiesSinglePhase.hpp b/opm/core/props/IncompPropertiesSinglePhase.hpp new file mode 100644 index 00000000..0be6df08 --- /dev/null +++ b/opm/core/props/IncompPropertiesSinglePhase.hpp @@ -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 . +*/ + +#ifndef OPM_INCOMPPROPERTIESSINGLEPHASE_HEADER_INCLUDED +#define OPM_INCOMPPROPERTIESSINGLEPHASE_HEADER_INCLUDED + + + +#include +#include +#include +#include + +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 diff --git a/opm/core/tof/FlowDiagnostics.cpp b/opm/core/tof/FlowDiagnostics.cpp new file mode 100644 index 00000000..a31ce498 --- /dev/null +++ b/opm/core/tof/FlowDiagnostics.cpp @@ -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 . +*/ + +#include + +#include +#include +#include + +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> computeFandPhi(const std::vector& pv, + const std::vector& ftof, + const std::vector& 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 D2; + std::vector 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 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 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& flowcap, + const std::vector& 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> computeSweep(const std::vector& flowcap, + const std::vector& 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 Ev; + std::vector 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 diff --git a/opm/core/tof/FlowDiagnostics.hpp b/opm/core/tof/FlowDiagnostics.hpp new file mode 100644 index 00000000..96b51234 --- /dev/null +++ b/opm/core/tof/FlowDiagnostics.hpp @@ -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 . +*/ + +#ifndef OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED +#define OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED + + +#include +#include + + +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> + computeFandPhi(const std::vector& pv, + const std::vector& ftof, + const std::vector& 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& flowcap, + const std::vector& 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> + computeSweep(const std::vector& flowcap, + const std::vector& storagecap); + +} // namespace Opm + +#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED diff --git a/tests/test_flowdiagnostics.cpp b/tests/test_flowdiagnostics.cpp new file mode 100644 index 00000000..a3d700e3 --- /dev/null +++ b/tests/test_flowdiagnostics.cpp @@ -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 . +*/ + +#include + +#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 +#include + +const std::vector pv(16, 18750.0); + +const std::vector 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 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 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 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 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 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 wrong_length(7, 0.0); + + + +using namespace Opm; + + +template +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); +}