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

Merge pull request #81 from qilicun/wpolymer_schedule

Browse Source 
add new polymerInflow constructor, get polymer injection rate from schedule section
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-12 09:38:34 +01:00
commit a20f0291a1
12 changed files with 173 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/sim_poly2p_comp_reorder.cpp
View File

@ -303,7 +303,7 @@ try
if (wells.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells()));
polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells(), simtimer.currentStepNum()));
} else {
polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,

2
examples/sim_poly2p_incomp_reorder.cpp
View File

@ -308,7 +308,7 @@ try
if (wells.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells()));
polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells(), simtimer.currentStepNum()));
} else {
polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,

77
opm/polymer/PolymerInflow.cpp
View File

@ -22,8 +22,15 @@
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/wells.h>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellPolymerProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellInjectionProperties.hpp>
#include <map>
#include <unordered_map>
#include <memory>
namespace Opm
{
@ -105,6 +112,74 @@ namespace Opm
}
}
void
PolymerInflowFromDeck::setInflowValues(Opm::DeckConstPtr deck,
size_t currentStep)
{
Opm::DeckKeywordConstPtr keyword = deck->getKeyword("WPOLYMER");
Schedule schedule(deck);
for (size_t recordNr = 0; recordNr < keyword->size(); recordNr++) {
DeckRecordConstPtr record = keyword->getRecord(recordNr);
const std::string& wellNamesPattern = record->getItem("WELL")->getTrimmedString(0);
std::string wellName = record->getItem("WELL")->getTrimmedString(0);
std::vector<WellPtr> wells = schedule.getWells(wellNamesPattern);
for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter) {
WellPtr well = *wellIter;
WellInjectionProperties injection = well->getInjectionProperties(currentStep);
if (injection.injectorType == WellInjector::WATER) {
WellPolymerProperties polymer = well->getPolymerProperties(currentStep);
wellPolymerRate_.insert(std::make_pair(wellName, polymer.m_polymerConcentration));
} else {
OPM_THROW(std::logic_error, "For polymer injector you must have a water injector");
}
}
}
}
/// Constructor.
/// @param[in] deck Input deck expected to contain WPOLYMER.
PolymerInflowFromDeck::PolymerInflowFromDeck(Opm::DeckConstPtr deck,
const Wells& wells,
const int num_cells,
size_t currentStep)
: sparse_inflow_(num_cells)
{
if (!deck->hasKeyword("WPOLYMER")) {
OPM_MESSAGE("PolymerInflowFromDeck initialized without WPOLYMER in current epoch.");
return;
}
setInflowValues(deck, currentStep);
std::unordered_map<std::string, double>::const_iterator map_it;
// Extract concentrations and put into cell->concentration map.
std::map<int, double> perfcell_conc;
for (size_t i = 0; i < wellPolymerRate_.size(); ++i) {
// Only use well name and polymer concentration.
// That is, we ignore salt concentration and group
// names.
int wix = 0;
for (; wix < wells.number_of_wells; ++wix) {
map_it = wellPolymerRate_.find(wells.name[wix]);
if (map_it == wellPolymerRate_.end()) {
OPM_THROW(std::runtime_error, "Could not find a match for well from WPOLYMER.");
} else {
break;
}
}
for (int j = wells.well_connpos[wix]; j < wells.well_connpos[wix+1]; ++j) {
const int perf_cell = wells.well_cells[j];
perfcell_conc[perf_cell] = map_it->second;
}
}
// Build sparse vector from map.
std::map<int, double>::const_iterator it = perfcell_conc.begin();
for (; it != perfcell_conc.end(); ++it) {
sparse_inflow_.addElement(it->second, it->first);
}
}
void PolymerInflowFromDeck::getInflowValues(const double /*step_start*/,
const double /*step_end*/,

19
opm/polymer/PolymerInflow.hpp
View File

@ -22,7 +22,12 @@
#include <opm/core/utility/SparseVector.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellPolymerProperties.hpp>
#include <vector>
#include <string>
#include <unordered_map>
struct Wells;
@ -90,6 +95,16 @@ namespace Opm
const Wells& wells,
const int num_cells);
/// Constructor.
/// \param[in] deck Input deck expected to contain WPOLYMER.
/// \param[in] wells Wells structure.
/// \param[in] num_cells Number of cells in grid.
/// \param[in] currentStep Number of current simulation step.
PolymerInflowFromDeck(Opm::DeckConstPtr deck,
const Wells& wells,
const int num_cells,
size_t currentStep);
/// Get inflow concentrations for all cells.
/// \param[in] step_start Start of timestep.
/// \param[in] step_end End of timestep.
@ -100,6 +115,10 @@ namespace Opm
std::vector<double>& poly_inflow_c) const;
private:
SparseVector<double> sparse_inflow_;
std::unordered_map<std::string, double> wellPolymerRate_;
void setInflowValues(Opm::DeckConstPtr deck,
size_t currentStep);
};

7
opm/polymer/PolymerProperties.cpp
View File

@ -135,12 +135,7 @@ namespace Opm
if (ads_index_ == Desorption) {
simpleAdsorptionBoth(c, c_ads, dc_ads_dc, if_with_der);
} else if (ads_index_ == NoDesorption) {
if (c < cmax) {
simpleAdsorption(cmax, c_ads);
dc_ads_dc = 0.;
} else {
simpleAdsorptionBoth(c, c_ads, dc_ads_dc, if_with_der);
}
simpleAdsorptionBoth(std::max(c, cmax), c_ads, dc_ads_dc, if_with_der);
} else {
OPM_THROW(std::runtime_error, "Invalid Adsoption index");
}

5
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
View File

@ -257,8 +257,7 @@ namespace Opm {
const SolutionState& state);
void
computeCmax(PolymerBlackoilState& state,
const ADB& c);
computeCmax(PolymerBlackoilState& state);
ADB
computeMc(const SolutionState& state) const;
@ -351,7 +350,7 @@ namespace Opm {
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
bool getConvergence(const double dt);
bool getConvergence(const double dt, const int it);
void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol,

89
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
View File

@ -275,15 +275,16 @@ namespace {
{
const V pvdt = geo_.poreVolume() / dt;
// Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
if (active_[Gas]) { updatePrimalVariableFromState(x); }
{
const SolutionState state = constantState(x, xw);
computeCmax(x, state.concentration);
computeAccum(state, 0);
computeWellConnectionPressures(state, xw);
}
std::vector<std::vector<double>> residual_history;
assemble(pvdt, x, xw, polymer_inflow);
@ -291,16 +292,11 @@ namespace {
bool converged = false;
double omega = 1.;
const double r0 = residualNorm();
residual_history.push_back(residuals());
converged = getConvergence(dt);
int it = 0;
std::cout << "\nIteration Residual\n"
<< std::setw(9) << it << std::setprecision(9)
<< std::setw(18) << r0 << std::endl;
converged = getConvergence(dt, it);
const int sizeNonLinear = residual_.sizeNonLinear();
@ -327,21 +323,19 @@ namespace {
assemble(pvdt, x, xw, polymer_inflow);
const double r = residualNorm();
residual_history.push_back(residuals());
converged = getConvergence(dt);
++it;
it += 1;
std::cout << std::setw(9) << it << std::setprecision(9)
<< std::setw(18) << r << std::endl;
converged = getConvergence(dt, it);
}
if (!converged) {
std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
// OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
}
// Update max concentration.
computeCmax(x);
}
@ -647,14 +641,14 @@ namespace {
template<class T>
void FullyImplicitBlackoilPolymerSolver<T>::computeCmax(PolymerBlackoilState& state,
const ADB& c)
void FullyImplicitBlackoilPolymerSolver<T>::computeCmax(PolymerBlackoilState& state)
{
const int nc = AutoDiffGrid::numCells(grid_);
V tmp = V::Zero(nc);
for (int i = 0; i < nc; ++i) {
cmax_(i) = std::max(cmax_(i), c.value()(i));
tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
}
std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin());
std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
}
@ -1901,7 +1895,7 @@ namespace {
template<class T>
bool
FullyImplicitBlackoilPolymerSolver<T>::getConvergence(const double dt)
FullyImplicitBlackoilPolymerSolver<T>::getConvergence(const double dt, const int iteration)
{
const double tol_mb = 1.0e-7;
const double tol_cnv = 1.0e-3;
@ -1917,10 +1911,12 @@ namespace {
double CNVW = 0.;
double CNVO = 0.;
double CNVG = 0.;
double CNVP = 0.;
double RW_sum = 0.;
double RO_sum = 0.;
double RG_sum = 0.;
double RP_sum = 0.;
double BW_avg = 0.;
double BO_avg = 0.;
@ -1930,7 +1926,7 @@ namespace {
const int pos = pu.phase_pos[Water];
const ADB& tempBW = rq_[pos].b;
V BW = 1./tempBW.value();
V RW = residual_.material_balance_eq[Water].value();
V RW = residual_.material_balance_eq[pos].value();
BW_avg = BW.sum()/nc;
const V tempV = RW.abs()/pv;
@ -1943,7 +1939,7 @@ namespace {
const int pos = pu.phase_pos[Oil];
const ADB& tempBO = rq_[pos].b;
V BO = 1./tempBO.value();
V RO = residual_.material_balance_eq[Oil].value();
V RO = residual_.material_balance_eq[pos].value();
BO_avg = BO.sum()/nc;
const V tempV = RO.abs()/pv;
@ -1956,7 +1952,7 @@ namespace {
const int pos = pu.phase_pos[Gas];
const ADB& tempBG = rq_[pos].b;
V BG = 1./tempBG.value();
V RG = residual_.material_balance_eq[Gas].value();
V RG = residual_.material_balance_eq[pos].value();
BG_avg = BG.sum()/nc;
const V tempV = RG.abs()/pv;
@ -1964,13 +1960,28 @@ namespace {
RG_sum = RG.sum();
}
double tempValue = tol_mb * pvSum /dt;
if (has_polymer_) {
const ADB& tempBW = rq_[poly_pos_].b;
V BW = 1./tempBW.value();
V RP = residual_.material_balance_eq[poly_pos_].value();
BW_avg = BW.sum()/nc;
const V tempV = RP.abs()/pv;
bool converged_MB = (fabs(BW_avg*RW_sum) < tempValue)
&& (fabs(BO_avg*RO_sum) < tempValue)
&& (fabs(BG_avg*RG_sum) < tempValue);
CNVP = BW_avg * dt * tempV.maxCoeff();
RP_sum = RP.sum();
}
bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv);
const double mass_balance_residual_water = fabs(BW_avg*RW_sum) * dt / pvSum;
const double mass_balance_residual_oil = fabs(BO_avg*RO_sum) * dt / pvSum;
const double mass_balance_residual_gas = fabs(BG_avg*RG_sum) * dt / pvSum;
const double mass_balance_residual_polymer = fabs(BW_avg*RP_sum) * dt / pvSum;
bool converged_MB = (mass_balance_residual_water < tol_mb)
&& (mass_balance_residual_oil< tol_mb)
&& (mass_balance_residual_gas < tol_mb)
&& (mass_balance_residual_polymer < tol_mb);
bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv) && (CNVP < tol_cnv);
double residualWellFlux = residual_.well_flux_eq.value().matrix().template lpNorm<Eigen::Infinity>();
double residualWell = residual_.well_eq.value().matrix().template lpNorm<Eigen::Infinity>();
@ -1979,11 +1990,25 @@ namespace {
bool converged = converged_MB && converged_CNV && converged_Well;
#ifdef OPM_VERBOSE
std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl;
std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV
<< " converged_Well " << converged_Well << " converged " << converged << std::endl;
#endif
if (iteration == 0) {
std::cout << "\nIter MB(OIL) MB(WATER) MB(GAS) MB(POLYMER) CNVW CNVO CNVG CNVP WELL-FLOW WELL-CNTRL\n";
}
const std::streamsize oprec = std::cout.precision(3);
const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
std::cout << std::setw(4) << iteration
<< std::setw(11) << mass_balance_residual_water
<< std::setw(11) << mass_balance_residual_oil
<< std::setw(11) << mass_balance_residual_gas
<< std::setw(11) << mass_balance_residual_polymer
<< std::setw(11) << CNVW
<< std::setw(11) << CNVO
<< std::setw(11) << CNVG
<< std::setw(11) << CNVP
<< std::setw(11) << residualWellFlux
<< std::setw(11) << residualWell
<< std::endl;
std::cout.precision(oprec);
std::cout.flags(oflags);
return converged;
}

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

@ -202,9 +202,10 @@ namespace {
WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow)
{
// Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
const SolutionState state = constantState(x, xw);
computeCmax(x, state.concentration);
computeAccum(state, 0);
const double atol = 1.0e-12;
@ -241,6 +242,9 @@ namespace {
std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
// OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
}
// Update max concentration.
computeCmax(x);
}
@ -483,15 +487,15 @@ namespace {
void
FullyImplicitCompressiblePolymerSolver::
computeCmax(PolymerBlackoilState& state,
const ADB& c)
computeCmax(PolymerBlackoilState& state)
{
const int nc = grid_.number_of_cells;
for (int i = 0; i < nc; ++i) {
cmax_(i) = std::max(cmax_(i), c.value()(i));
V tmp = V::Zero(nc);
for (int i = 0; i < nc; ++i) {
tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
}
std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin());
std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
}

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

@ -196,8 +196,7 @@ namespace Opm {
const ADB& c) const;
void
computeCmax(PolymerBlackoilState& state,
const ADB& c);
computeCmax(PolymerBlackoilState& state);
ADB
computeMc(const SolutionState& state) const;

16
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
View File

@ -233,9 +233,8 @@ namespace Opm {
double res_factor = polymer_props_.resFactor();
double factor = (res_factor -1.) / max_ads;
V rk = one + factor * ads;
V krw_eff = krw / rk;
return krw_eff;
return krw / rk;
}
@ -257,17 +256,8 @@ namespace Opm {
double res_factor = polymer_props_.resFactor();
double factor = (res_factor - 1.) / max_ads;
ADB rk = one + ads * factor;
ADB dkrw_ds = krw / rk;
ADB dkrw_dc = -factor * krw / (rk * rk) * ads ;
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dkrw_ds.derivative()[block] * sw.derivative()[block]
+ dkrw_dc.derivative()[block] * c.derivative()[block];
}
return ADB::function(krw_eff, jacs);
return krw / rk;
}
}// namespace Opm

2
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
View File

@ -335,7 +335,7 @@ namespace Opm
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_)));
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,

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

@ -285,7 +285,7 @@ namespace Opm
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_)));
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,