OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added well management to SimulatorCompressiblePolymer.

Browse Source 
Another polymer-specific overload of computeFractionalFlow()
has been added in support of this (this overload deals with
compressible properties.

Also added pressure normalization for situations with arbitrary
absolute pressure.
This commit is contained in:
Atgeirr Flø Rasmussen
2012-10-15 14:23:13 +02:00
parent f5861f0cc8
commit 2b292108ad
9 changed files with 172 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
examples/sim_poly2p_comp_reorder.cpp
View File

@@ -76,7 +76,7 @@ main(int argc, char** argv)
{
using namespace Opm;
std::cout << "\n================ Test program for incompressible two-phase flow with polymer ===============\n\n";
std::cout << "\n================ Test program for weakly compressible two-phase flow with polymer ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
@@ -241,12 +241,13 @@ main(int argc, char** argv)
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
param.getDefault("poly_amount", poly_props.cMax()));
WellsManager wells;
SimulatorCompressiblePolymer simulator(param,
*grid->c_grid(),
*props,
poly_props,
rock_comp->isActive() ? rock_comp.get() : 0,
0, // wells
wells,
polymer_inflow,
src,
bcs.c_bcs(),
@@ -324,7 +325,7 @@ main(int argc, char** argv)
*props,
poly_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells.c_wells(),
wells,
*polymer_inflow,
src,
bcs.c_bcs(),

98
opm/polymer/SimulatorCompressiblePolymer.cpp
View File

@@ -17,7 +17,6 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
@@ -39,6 +38,7 @@
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
@@ -52,7 +52,7 @@
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/polymerUtilities.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
@@ -92,7 +92,7 @@ namespace Opm
const BlackoilPropertiesInterface& props,
const PolymerProperties& poly_props,
const RockCompressibility* rock_comp_props,
const Wells* wells,
WellsManager& wells_manager,
const PolymerInflowInterface& polymer_inflow,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
@@ -111,6 +111,9 @@ namespace Opm
bool output_vtk_;
std::string output_dir_;
int output_interval_;
// Parameters for well control
bool check_well_controls_;
int max_well_control_iterations_;
// Parameters for transport solver.
int num_transport_substeps_;
bool use_segregation_split_;
@@ -119,11 +122,11 @@ namespace Opm
const BlackoilPropertiesInterface& props_;
const PolymerProperties& poly_props_;
const RockCompressibility* rock_comp_props_;
WellsManager& wells_manager_;
const Wells* wells_;
const PolymerInflowInterface& polymer_inflow_;
const std::vector<double>& src_;
const FlowBoundaryConditions* bcs_;
const LinearSolverInterface& linsolver_;
const double* gravity_;
// Solvers
CompressibleTpfaPolymer psolver_;
@@ -142,7 +145,7 @@ namespace Opm
const BlackoilPropertiesInterface& props,
const PolymerProperties& poly_props,
const RockCompressibility* rock_comp_props,
const Wells* wells,
WellsManager& wells_manager,
const PolymerInflowInterface& polymer_inflow,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
@@ -150,7 +153,7 @@ namespace Opm
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props,
wells, polymer_inflow, src, bcs, linsolver, gravity));
wells_manager, polymer_inflow, src, bcs, linsolver, gravity));
}
@@ -173,7 +176,7 @@ namespace Opm
const BlackoilPropertiesInterface& props,
const PolymerProperties& poly_props,
const RockCompressibility* rock_comp_props,
const Wells* wells,
WellsManager& wells_manager,
const PolymerInflowInterface& polymer_inflow,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
@@ -183,18 +186,18 @@ namespace Opm
props_(props),
poly_props_(poly_props),
rock_comp_props_(rock_comp_props),
wells_(wells),
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
polymer_inflow_(polymer_inflow),
src_(src),
bcs_(bcs),
linsolver_(linsolver),
gravity_(gravity),
psolver_(grid, props, rock_comp_props, poly_props, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells),
tsolver_(grid, props, poly_props, *rock_comp_props,
gravity, wells_manager.c_wells()),
tsolver_(grid, props, poly_props,
TransportModelCompressiblePolymer::Bracketing,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30))
@@ -215,6 +218,10 @@ namespace Opm
output_interval_ = param.getDefault("output_interval", 1);
}
// Well control related init.
check_well_controls_ = param.getDefault("check_well_controls", false);
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Transport related init.
TransportModelCompressiblePolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
@@ -228,7 +235,7 @@ namespace Opm
tsolver_.setPreferredMethod(method);
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity_ != 0 && use_segregation_split_){
if (gravity != 0 && use_segregation_split_) {
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
@@ -261,7 +268,6 @@ namespace Opm
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
@@ -301,15 +307,75 @@ namespace Opm
initial_pressure = state.pressure();
// Solve pressure.
// Solve pressure equation.
if (check_well_controls_) {
computeFractionalFlow(props_, poly_props_, allcells_,
state.pressure(), state.surfacevol(), state.saturation(),
state.concentration(), state.maxconcentration(),
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
bool well_control_passed = !check_well_controls_;
int well_control_iteration = 0;
do {
// Run solver
pressure_timer.start();
psolver_.solve(timer.currentStepLength(), state, well_state);
// Renormalize pressure if both fluids and rock are
// incompressible, and there are no pressure
// conditions (bcs or wells). It is deemed sufficient
// for now to renormalize using geometric volume
// instead of pore volume.
if (psolver_.singularPressure()) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.0;
double av_press = 0.0;
double tot_vol = 0.0;
const int num_cells = grid_.number_of_cells;
for (int cell = 0; cell < num_cells; ++cell) {
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
tot_vol += grid_.cell_volumes[cell];
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
}
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
for (int well = 0; well < num_wells; ++well) {
well_state.bhp()[well] += ren_const;
}
}
// Stop timer and report
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
} while (false);
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -358,7 +424,7 @@ namespace Opm
polyprod += substep_polyprod;
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, stepsize,
state.saturation(), state.surfacevol(),
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());
}
}

23
opm/polymer/SimulatorCompressiblePolymer.hpp
View File

@@ -33,6 +33,7 @@ namespace Opm
class BlackoilPropertiesInterface;
class PolymerProperties;
class RockCompressibility;
class WellsManager;
class PolymerInflowInterface;
class LinearSolverInterface;
class SimulatorTimer;
@@ -60,22 +61,22 @@ namespace Opm
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] poly_props polymer properties
/// \param[in] rock_comp if non-null, rock compressibility properties
/// \param[in] wells if non-null, wells data structure
/// \param[in] polymer_inflow polymer inflow controls
/// \param[in] src source terms
/// \param[in] bcs boundary conditions, treat as all noflow if null
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] poly_props polymer properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells_manager well manager, may manage no (null) wells
/// \param[in] polymer_inflow polymer inflow controls
/// \param[in] src source terms
/// \param[in] bcs boundary conditions, treat as all noflow if null
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorCompressiblePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const PolymerProperties& poly_props,
const RockCompressibility* rock_comp_props,
const Wells* wells,
WellsManager& wells_manager,
const PolymerInflowInterface& polymer_inflow,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,

6
opm/polymer/SimulatorPolymer.cpp
View File

@@ -52,7 +52,7 @@
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/polymerUtilities.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
@@ -129,7 +129,6 @@ namespace Opm
const PolymerInflowInterface& polymer_inflow_;
const std::vector<double>& src_;
const FlowBoundaryConditions* bcs_;
const LinearSolverInterface& linsolver_;
const double* gravity_;
// Solvers
IncompTpfaPolymer psolver_;
@@ -194,7 +193,6 @@ namespace Opm
polymer_inflow_(polymer_inflow),
src_(src),
bcs_(bcs),
linsolver_(linsolver),
gravity_(gravity),
psolver_(grid, props, rock_comp_props, poly_props, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
@@ -238,7 +236,7 @@ namespace Opm
tsolver_.setPreferredMethod(method);
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
if (gravity != 0 && use_segregation_split_) {
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}

21
opm/polymer/SimulatorPolymer.hpp
View File

@@ -61,17 +61,16 @@ namespace Opm
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] poly_props polymer properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells if non-null, wells data structure
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] polymer_inflow polymer inflow controls
/// \param[in] src source terms
/// \param[in] bcs boundary conditions, treat as all noflow if null
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] poly_props polymer properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells_manager well manager, may manage no (null) wells
/// \param[in] polymer_inflow polymer inflow controls
/// \param[in] src source terms
/// \param[in] bcs boundary conditions, treat as all noflow if null
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorPolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,

2
opm/polymer/TransportModelCompressiblePolymer.cpp
View File

@@ -154,14 +154,12 @@ namespace Opm
TransportModelCompressiblePolymer::TransportModelCompressiblePolymer(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const PolymerProperties& polyprops,
const RockCompressibility& rock_comp,
const SingleCellMethod method,
const double tol,
const int maxit)
: grid_(grid),
props_(props),
polyprops_(polyprops),
rock_comp_(rock_comp),
darcyflux_(0),
porevolume0_(0),
porevolume_(0),

3
opm/polymer/TransportModelCompressiblePolymer.hpp
View File

@@ -20,7 +20,6 @@
#ifndef OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
@@ -63,7 +62,6 @@ namespace Opm
TransportModelCompressiblePolymer(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const PolymerProperties& polyprops,
const RockCompressibility& rock_comp,
const SingleCellMethod method,
const double tol,
const int maxit);
@@ -134,7 +132,6 @@ namespace Opm
const UnstructuredGrid& grid_;
const BlackoilPropertiesInterface& props_;
const PolymerProperties& polyprops_;
const RockCompressibility& rock_comp_;
const double* darcyflux_; // one flux per grid face
const double* porevolume0_; // one volume per cell
const double* porevolume_; // one volume per cell

42
opm/polymer/polymerUtilities.cpp
View File

@@ -127,6 +127,48 @@ namespace Opm
}
/// Computes the fractional flow for each cell in the cells argument
/// @param[in] props rock and fluid properties
/// @param[in] polyprops polymer properties
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases)
/// @param[in] c concentration values
/// @param[in] cmax max polymer concentration experienced by cell
/// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const Opm::PolymerProperties& polyprops,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& c,
const std::vector<double>& cmax,
std::vector<double>& fractional_flows)
{
int num_cells = cells.size();
int num_phases = props.numPhases();
if (num_phases != 2) {
THROW("computeFractionalFlow() assumes 2 phases.");
}
fractional_flows.resize(num_cells*num_phases);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
std::vector<double> mu(num_cells*num_phases);
props.viscosity(num_phases, &p[0], &z[0], &cells[0], &mu[0], 0);
double mob[2]; // here we assume num_phases=2
for (int cell = 0; cell < num_cells; ++cell) {
double* kr_cell = &kr[2*cell];
double* mu_cell = &mu[2*cell];
polyprops.effectiveMobilities(c[cell], cmax[cell], mu_cell, kr_cell, mob);
fractional_flows[2*cell] = mob[0] / (mob[0] + mob[1]);
fractional_flows[2*cell + 1] = mob[1] / (mob[0] + mob[1]);
}
}
/// @brief Computes injected and produced volumes of all phases,
/// and injected and produced polymer mass.
/// Note 1: assumes that only the first phase is injected.

20
opm/polymer/polymerUtilities.hpp
View File

@@ -85,6 +85,26 @@ namespace Opm
const std::vector<double>& cmax,
std::vector<double>& fractional_flows);
/// Computes the fractional flow for each cell in the cells argument
/// @param[in] props rock and fluid properties
/// @param[in] polyprops polymer properties
/// @param[in] cells cells with which the saturation values are associated
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all phases)
/// @param[in] c concentration values
/// @param[in] cmax max polymer concentration experienced by cell
/// @param[out] fractional_flow the fractional flow for each phase for each cell.
void computeFractionalFlow(const Opm::BlackoilPropertiesInterface& props,
const Opm::PolymerProperties& polyprops,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& c,
const std::vector<double>& cmax,
std::vector<double>& fractional_flows);
/// @brief Computes injected and produced volumes of all phases,
/// and injected and produced polymer mass.
/// Note 1: assumes that only the first phase is injected.