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

Merge pull request #752 from atgeirr/compute-well-pairs

Browse Source 
Improvements for flow diagnostic tools
This commit is contained in:
Atgeirr Flø Rasmussen 2015-02-15 12:49:29 +01:00
commit 1315be7be6
3 changed files with 157 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
examples/compute_tof.cpp
View File

@ -39,6 +39,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/pressure/IncompTpfaSinglePhase.hpp>
#include <opm/core/flowdiagnostics/FlowDiagnostics.hpp>
#include <opm/core/flowdiagnostics/TofReorder.hpp>
#include <opm/core/flowdiagnostics/TofDiscGalReorder.hpp>
@ -203,14 +204,6 @@ try
}
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.
bool output = param.getDefault("output", true);
std::ofstream epoch_os;
@ -252,68 +245,96 @@ try
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<double> src(num_cells, 0.0);
std::vector<double> transport_src;
computeTransportSourceSinglePhase(grid, src, flux, 1.0,
&wells, wellrates, transport_src);
// Solve time-of-flight.
transport_timer.start();
std::vector<double> tof;
std::vector<double> tracer;
Opm::SparseTable<int> tracerheads;
if (compute_tracer) {
buildTracerheadsFromWells(wells, start_from_injectors, tracerheads);
}
if (use_dg) {
if (compute_tracer) {
dg_solver->solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else {
dg_solver->solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
}
} else {
Opm::TofReorder tofsolver(grid, use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else {
tofsolver.solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
total_timer.stop();
std::string tof_filenames[2] = { output_dir + "/ftof.txt", output_dir + "/btof.txt" };
std::string tracer_filenames[2] = { output_dir + "/ftracer.txt", output_dir + "/btracer.txt" };
std::vector<double> tracers[2];
// Output.
if (output) {
std::string tof_filename = output_dir + "/tof.txt";
std::ofstream tof_stream(tof_filename.c_str());
tof_stream.precision(16);
std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
// We compute tof twice, direction == 0 is from injectors, 1 is from producers.
for (int direction = 0; direction < 2; ++direction) {
// Turn direction of flux and flip source terms if starting from producers.
if (direction == 1) {
for (auto it = flux.begin(); it != flux.end(); ++it) {
(*it) = -(*it);
}
for (auto it = transport_src.begin(); it != transport_src.end(); ++it) {
(*it) = -(*it);
}
}
// Solve time-of-flight.
transport_timer.start();
std::vector<double> tof;
std::vector<double> tracer;
Opm::SparseTable<int> tracerheads;
if (compute_tracer) {
std::string tracer_filename = output_dir + "/tracer.txt";
std::ofstream tracer_stream(tracer_filename.c_str());
tracer_stream.precision(16);
const int nt = tracer.size()/num_cells;
for (int i = 0; i < nt*num_cells; ++i) {
tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
buildTracerheadsFromWells(wells, direction == 0, tracerheads);
}
if (use_dg) {
if (compute_tracer) {
dg_solver->solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else {
dg_solver->solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
}
} else {
Opm::TofReorder tofsolver(grid, use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(flux.data(), porevol.data(), transport_src.data(), tracerheads, tof, tracer);
} else {
tofsolver.solveTof(flux.data(), porevol.data(), transport_src.data(), tof);
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
if (direction == 0) {
std::cout << "Forward ";
} else {
std::cout << "Backward ";
}
std::cout << "time-of-flight/tracer solve took: " << tt << " seconds." << std::endl;
ttime += tt;
// Output.
if (output) {
std::string tof_filename = tof_filenames[direction];
std::ofstream tof_stream(tof_filename.c_str());
tof_stream.precision(16);
std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
if (compute_tracer) {
std::string tracer_filename = tracer_filenames[direction];
std::ofstream tracer_stream(tracer_filename.c_str());
tracer_stream.precision(16);
const int nt = tracer.size()/num_cells;
for (int i = 0; i < nt*num_cells; ++i) {
tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
}
tracers[direction] = tracer;
}
}
}
// If we have tracers, compute well pairs.
if (compute_tracer) {
auto wp = Opm::computeWellPairs(wells, porevol, tracers[0], tracers[1]);
std::string wellpair_filename = output_dir + "/wellpairs.txt";
std::ofstream wellpair_stream(wellpair_filename.c_str());
const int nwp = wp.size();
for (int ii = 0; ii < nwp; ++ii) {
wellpair_stream << std::get<0>(wp[ii]) << ' ' << std::get<1>(wp[ii]) << ' ' << std::get<2>(wp[ii]) << '\n';
}
}
total_timer.stop();
std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Tof/tracer time: " << ttime << std::endl;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";

60
opm/core/flowdiagnostics/FlowDiagnostics.cpp
View File

@ -18,6 +18,7 @@
*/
#include <opm/core/flowdiagnostics/FlowDiagnostics.hpp>
#include <opm/core/wells.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
@ -163,4 +164,63 @@ namespace Opm
/// \brief Compute volumes associated with injector-producer pairs.
///
/// \param[in] wells wells structure, containing NI injector wells and NP producer wells.
/// \param[in] porevol pore volume of each grid cell
/// \param[in] ftracer array of forward (injector) tracer values, NI per cell
/// \param[in] btracer array of backward (producer) tracer values, NP per cell
/// \return a vector of tuples, one tuple for each injector-producer pair,
/// where the first and second elements are well indices for the
/// injector and producer, and the third element is the pore volume
/// associated with that pair.
std::vector<std::tuple<int, int, double> >
computeWellPairs(const Wells& wells,
const std::vector<double>& porevol,
const std::vector<double>& ftracer,
const std::vector<double>& btracer)
{
// Identify injectors and producers.
std::vector<int> inj;
std::vector<int> prod;
const int nw = wells.number_of_wells;
for (int w = 0; w < nw; ++w) {
if (wells.type[w] == INJECTOR) {
inj.push_back(w);
} else {
prod.push_back(w);
}
}
// Check sizes of input arrays.
const int nc = porevol.size();
if (nc * inj.size() != ftracer.size()) {
OPM_THROW(std::runtime_error, "computeWellPairs(): wrong size of input array ftracer.");
}
if (nc * prod.size() != btracer.size()) {
OPM_THROW(std::runtime_error, "computeWellPairs(): wrong size of input array btracer.");
}
// Compute associated pore volumes.
std::vector<std::tuple<int, int, double> > result;
const int num_inj = inj.size();
const int num_prod = prod.size();
for (int inj_ix = 0; inj_ix < num_inj; ++inj_ix) {
for (int prod_ix = 0; prod_ix < num_prod; ++prod_ix) {
double assoc_porevol = 0.0;
for (int c = 0; c < nc; ++c) {
assoc_porevol += porevol[c]
* ftracer[num_inj * c + inj_ix]
* btracer[num_prod * c + prod_ix];
}
result.push_back(std::make_tuple(inj[inj_ix], prod[prod_ix], assoc_porevol));
}
}
return result;
}
} // namespace Opm

19
opm/core/flowdiagnostics/FlowDiagnostics.hpp
View File

@ -23,7 +23,9 @@
#include <vector>
#include <utility>
#include <tuple>
struct Wells;
namespace Opm
{
@ -79,6 +81,23 @@ namespace Opm
computeSweep(const std::vector<double>& flowcap,
const std::vector<double>& storagecap);
/// \brief Compute volumes associated with injector-producer pairs.
///
/// \param[in] wells wells structure, containing NI injector wells and NP producer wells.
/// \param[in] porevol pore volume of each grid cell
/// \param[in] ftracer array of forward (injector) tracer values, NI per cell
/// \param[in] btracer array of backward (producer) tracer values, NP per cell
/// \return a vector of tuples, one tuple for each injector-producer pair,
/// where the first and second elements are well indices for the
/// injector and producer, and the third element is the pore volume
/// associated with that pair.
std::vector<std::tuple<int, int, double>>
computeWellPairs(const Wells& wells,
const std::vector<double>& porevol,
const std::vector<double>& ftracer,
const std::vector<double>& btracer);
} // namespace Opm
#endif // OPM_FLOWDIAGNOSTICS_HEADER_INCLUDED