OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-20 01:32:57 -06:00
Code Issues Packages Projects Releases Wiki Activity

output watercut by using utilities.cpp/computeInjectedProduced function.

Browse Source
This commit is contained in:
Liu Ming 2014-01-22 14:17:04 +08:00
parent 8987a4c1d7
commit 12318fe754
6 changed files with 334 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -142,7 +142,6 @@ try
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
// Linear solver.
LinearSolverFactory linsolver(param);

27
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -163,7 +163,8 @@ namespace {
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver)
const LinearSolverInterface& linsolver
)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
@ -192,7 +193,8 @@ namespace {
step(const double dt,
PolymerBlackoilState& x ,
WellState& xw,
const std::vector<double>& polymer_inflow)
const std::vector<double>& polymer_inflow,
std::vector<double>& src)
{
const SolutionState state = constantState(x, xw);
@ -202,8 +204,7 @@ namespace {
const double atol = 1.0e-12;
const double rtol = 5.0e-8;
const int maxit = 15;
assemble(dt, x, xw, polymer_inflow);
assemble(dt, x, xw, polymer_inflow, src);
const double r0 = residualNorm();
const double r_polymer = residual_.mass_balance[2].value().matrix().lpNorm<Eigen::Infinity>();
@ -217,7 +218,7 @@ namespace {
const V dx = solveJacobianSystem();
updateState(dx, x, xw);
assemble(dt, x, xw, polymer_inflow);
assemble(dt, x, xw, polymer_inflow, src);
const double r = residualNorm();
@ -482,9 +483,11 @@ namespace {
assemble(const double dt,
const PolymerBlackoilState& x ,
const WellState& xw,
const std::vector<double>& polymer_inflow)
const std::vector<double>& polymer_inflow,
std::vector<double>& src)
{
// Create the primary variables.
//
const SolutionState state = variableState(x, xw);
const V pvdt = geo_.poreVolume() / dt;
// -------- Mass balance equations --------
@ -519,7 +522,6 @@ namespace {
// Add the extra (flux) terms to the gas mass balance equations
// from gas dissolved in the oil phase.
// The extra terms in the accumulation part of the equation are already handled.
// -------- Well equation, and well contributions to the mass balance equations --------
// Contribution to mass balance will have to wait.
@ -528,6 +530,9 @@ namespace {
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
for (int i = 0; i < nc; ++i) {
src[i] = 0.0;
}
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
const V transw = Eigen::Map<const V>(wells_.WI, nperf);
@ -603,6 +608,9 @@ namespace {
well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
// DUMP(well_contribs[phase]);
residual_.mass_balance[phase] += well_contribs[phase];
for (int cell = 0; cell < nc; ++cell) {
src[cell] += well_contribs[phase].value()[cell];
}
}
// well rates contribs to polymer mass balance eqn.
@ -647,6 +655,11 @@ namespace {
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// for (int i = 0; i < nc; ++i) {
// std::cout << src[i] << " ";
// if ((i+1) % 10 == 0)
// std::cout<<std::endl;
// }
// DUMP(residual_.well_eq);
}

6
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -79,7 +79,8 @@ namespace Opm {
step(const double dt ,
PolymerBlackoilState& state ,
WellState& wstate,
const std::vector<double>& polymer_inflow);
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
private:
// Types and enums
@ -160,7 +161,8 @@ namespace Opm {
assemble(const double dt,
const PolymerBlackoilState& x,
const WellState& xw,
const std::vector<double>& polymer_inflow);
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
V solveJacobianSystem() const;

75
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

@ -29,7 +29,7 @@
#include <opm/polymer/fullyimplicit/GeoProps.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
#include <opm/polymer/fullyimplicit/utilities.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
@ -162,6 +162,7 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@ -178,6 +179,7 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
dm["surfvolume"] = &state.surfacevol();
std::vector<double> cell_velocity;
@ -205,6 +207,10 @@ namespace Opm
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
#if 0
//well rate should be modified by using effective viscosity
//and effective relperm for water
static void outputWellStateMatlab(const Opm::WellState& well_state,
const int step,
const std::string& output_dir)
@ -234,9 +240,9 @@ namespace Opm
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
#endif
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
@ -244,7 +250,7 @@ namespace Opm
if (!os) {
OPM_THROW(std::runtime_error, "Failed to open " << fname);
}
watercut.write(os);
watercut.write(os);
}
#if 0
@ -331,25 +337,26 @@ namespace Opm
} else {
computePorevolume(grid_, props_.porosity(), porevol);
}
// const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
std::vector<double> polymer_inflow_c(grid_.number_of_cells);
std::vector<double> transport_src(grid_.number_of_cells);
// Main simulation loop.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
#if 0
// These must be changed for three-phase.
double init_surfvol[2] = { 0.0 };
double inplace_surfvol[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
#endif
std::vector<double> fractional_flows;
@ -376,7 +383,7 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
// outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
}
@ -392,14 +399,15 @@ namespace Opm
bool well_control_passed = !check_well_controls_;
int well_control_iteration = 0;
do {
// Process transport sources (to include bdy terms and well flows).
// Opm::computeTransportSource(props_, wells_, well_state, transport_src);
// Run solver.
const double current_time = timer.currentTime();
double stepsize = timer.currentStepLength();
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
solver_timer.start();
std::vector<double> initial_pressure = state.pressure();
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
// Stop timer and report.
solver_timer.stop();
const double st = solver_timer.secsSinceStart();
@ -428,13 +436,44 @@ namespace Opm
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
if (rock_comp_props_ && rock_comp_props_->isActive()) {
initial_porevol = porevol;
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double polyinj = 0;
double polyprod = 0;
Opm::computeInjectedProduced(props_, polymer_props_,
state,
transport_src, polymer_inflow_c, timer.currentStepLength(),
injected, produced,
polyinj, polyprod);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
watercut.push(timer.currentTime() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
std::cout.precision(5);
const int width = 18;
std::cout << "\nMass balance report.\n";
std::cout << " Injected reservoir volumes: "
<< std::setw(width) << injected[0]
<< std::setw(width) << injected[1] << std::endl;
std::cout << " Produced reservoir volumes: "
<< std::setw(width) << produced[0]
<< std::setw(width) << produced[1] << std::endl;
std::cout << " Total inj reservoir volumes: "
<< std::setw(width) << tot_injected[0]
<< std::setw(width) << tot_injected[1] << std::endl;
std::cout << " Total prod reservoir volumes: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1] << std::endl;
// The reports below are geared towards two phases only.
#if 0
// Report mass balances.
@ -443,10 +482,6 @@ namespace Opm
Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
injected, produced);
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
std::cout.precision(5);
const int width = 18;
std::cout << "\nMass balance report.\n";
@ -493,9 +528,9 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
#if 0
// outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
outputWaterCut(watercut, output_dir_);
#if 0
if (wells_) {
outputWellReport(wellreport, output_dir_);
}

174
opm/polymer/fullyimplicit/utilities.cpp Normal file
View File

@ -0,0 +1,174 @@
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/linalg/blas_lapack.h>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/polymer/fullyimplicit/utilities.hpp>
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <iostream>
#include <iterator>
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (reservoir volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropsAdInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src)
{
int nc = props.numCells();
transport_src.clear();
transport_src.resize(nc, 0.0);
// Well contributions.
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
if (np != 2) {
OPM_THROW(std::runtime_error, "computeTransportSource() requires a 2 phase case.");
}
std::vector<double> A(np*np);
for (int w = 0; w < nw; ++w) {
const double* comp_frac = wells->comp_frac + np*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_state.perfRates()[perf];
if (perf_rate > 0.0) {
// perf_rate is a total inflow reservoir rate, we want a surface water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Reservoir rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
perf_rate *= comp_frac[0]; // Water reservoir volume rate.
}
}
transport_src[perf_cell] += perf_rate;
}
}
}
}
/// @brief Computes injected and produced volumes of all phases,
/// and injected and produced polymer mass - in the compressible case.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] polyprops polymer properties
/// @param[in] state state variables (pressure, fluxes etc.)
/// @param[in] transport_src if < 0: total reservoir volume outflow,
/// if > 0: first phase *surface volume* inflow.
/// @param[in] inj_c injected concentration by cell
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/transport_src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// @param[out] polyinj injected mass of polymer
/// @param[out] polyprod produced mass of polymer
void computeInjectedProduced(const BlackoilPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerBlackoilState& state,
const std::vector<double>& transport_src,
const std::vector<double>& inj_c,
const double dt,
double* injected,
double* produced,
double& polyinj,
double& polyprod)
{
const int num_cells = transport_src.size();
if (props.numCells() != num_cells) {
OPM_THROW(std::runtime_error, "Size of transport_src vector does not match number of cells in props.");
}
const int np = props.numPhases();
if (int(state.saturation().size()) != num_cells*np) {
OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
}
std::vector<int> cells(num_cells);
const V p = Eigen::Map<const V>(&state.pressure()[0], num_cells, 1);
const DataBlock s = Eigen::Map<const DataBlock>(&state.saturation()[0], num_cells, np);
const V sw = s.col(0);
const V so = s.col(1);
const V c = Eigen::Map<const V>(&state.concentration()[0], num_cells, 1);
const V cmax = Eigen::Map<const V>(&state.maxconcentration()[0], num_cells, 1);
const V trans_src = Eigen::Map<const V>(&transport_src[0], num_cells, 1);
V src = V::Constant(num_cells, -1.0); // negative is injec, positive is producer.
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
if(transport_src[cell] > 0.0) {
src[cell] = 1.0;
}
}
const Selector<double> src_selector(src);
const V one = V::Constant(num_cells, 1.0);
const V zero = V::Zero(num_cells);
const std::vector<V> kr = props.relperm(sw, so, zero, cells);
const V muw = props.muWat(p, cells);
const V muo = props.muOil(p, zero, cells);
const V krw_eff = polymer_props.effectiveRelPerm(c, cmax, kr[0]);
const V inv_muw_eff = polymer_props.effectiveInvWaterVisc(c, muw.data());
std::vector<V> mob(np);
mob[0] = krw_eff * inv_muw_eff;
mob[1] = kr[1] / muo;
const V watmob_c = src_selector.select(mob[0], one);
const V oilmob_c = src_selector.select(mob[1], zero);
const V flux = trans_src * dt;
const V totmob_c = watmob_c + oilmob_c;
const V wat_src = flux * (watmob_c / totmob_c);
const V oil_src = flux * (oilmob_c / totmob_c);
const V mc = polymer_props.polymerWaterVelocityRatio(c);
polyinj = 0.0;
polyprod = 0.0;
std::fill(injected, injected + np , 0.0);
std::fill(produced, produced + np , 0.0);
for (int cell = 0; cell < num_cells; ++cell) {
if (wat_src[cell] < 0) {
injected[0] += wat_src[cell];
polyinj += injected[0] * inj_c[cell];
} else {
produced[0] += wat_src[cell];
produced[1] += oil_src[cell];
polyprod += produced[0] * mc[cell];
}
}
}
}

81
opm/polymer/fullyimplicit/utilities.hpp Normal file
View File

@ -0,0 +1,81 @@
#ifndef OPM_UTILITIES_HEADER_INCLUDED
#define OPM_UTILITIES_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPropsAdInterface.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <algorithm>
#include <cmath>
#include <functional>
#include <limits>
#include <iostream>
#include <iterator>
namespace Opm
{
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (reservoir volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropsAdInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src);
/// @brief Computes injected and produced volumes of all phases,
/// and injected and produced polymer mass - in the compressible case.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] polyprops polymer properties
/// @param[in] state state variables (pressure, fluxes etc.)
/// @param[in] transport_src if < 0: total reservoir volume outflow,
/// if > 0: first phase *surface volume* inflow.
/// @param[in] inj_c injected concentration by cell
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/transport_src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// @param[out] polyinj injected mass of polymer
/// @param[out] polyprod produced mass of polymer
void computeInjectedProduced(const BlackoilPropsAdInterface& props,
const Opm::PolymerPropsAd& polymer_props,
const PolymerBlackoilState& state,
const std::vector<double>& transport_src,
const std::vector<double>& inj_c,
const double dt,
double* injected,
double* produced,
double& polyinj,
double& polyprod);
} //namespace Opm
#endif //OPM_UTILITIES_HEADER_INCLUDED