make some tricks for initialing polymer_inflow, this maybe wrong, but now no bugs.

This commit is contained in:
Liu Ming
2014-10-09 14:23:36 +08:00
parent 73c2032974
commit c2d6347271
5 changed files with 54 additions and 49 deletions

View File

@@ -75,7 +75,7 @@ namespace Opm
const std::vector<double>& facepressure() const { return state_blackoil_.facepressure(); }
const std::vector<double>& faceflux () const { return state_blackoil_.faceflux(); }
const std::vector<double>& saturation () const { return state_blackoil_.saturation(); }
const std::vector<double>& gasoilration() const { return state_blackoil_.gasoilratio(); }
const std::vector<double>& gasoilratio() const { return state_blackoil_.gasoilratio(); }
const std::vector<double>& rv () const { return state_blackoil_.rv(); }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& maxconcentration() const { return cmax_; }

View File

@@ -195,7 +195,7 @@ namespace {
, cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
, ops_ (grid)
, wops_ (wells)
, cmax_(V::Zero(numCells(grid)))
, cmax_(V::Zero(Opm::AutoDiffGrid::numCells(grid)))
, has_disgas_(has_disgas)
, has_vapoil_(has_vapoil)
, has_polymer_(has_polymer)
@@ -280,7 +280,7 @@ namespace {
std::vector<std::vector<double>> residual_history;
assemble(pvdt, x, xw);
assemble(pvdt, x, xw, polymer_inflow);
bool converged = false;
@@ -319,7 +319,7 @@ namespace {
updateState(dx, x, xw);
assemble(pvdt, x, xw);
assemble(pvdt, x, xw, polymer_inflow);
const double r = residualNorm();
@@ -621,7 +621,7 @@ namespace {
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const double rho_rock = polymer_props_ad_.rockDensity();
const V phi = Eigen::Map<const V>(& fluid_.porosity()[0], numCells(grid_), 1);
const V phi = Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1);
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
// compute total phases and determin polymer position.
rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads;
@@ -636,7 +636,7 @@ namespace {
void FullyImplicitBlackoilPolymerSolver<T>::computeCmax(PolymerBlackoilState& state,
const ADB& c)
{
const int nc = numCells(grid_);
const int nc = AutoDiffGrid::numCells(grid_);
for (int i = 0; i < nc; ++i) {
cmax_(i) = std::max(cmax_(i), c.value()(i));
}
@@ -795,7 +795,6 @@ namespace {
}
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
// Add polymer equation.
if (has_polymer_) {
residual_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
@@ -1339,10 +1338,10 @@ namespace {
const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
varstart += dxvar.size();
const V dc = null;
if (has_polymer_) {
const V dc = subset(dx, Span(nc, 1, varstart));
} else {
const V dc = null;
varstart += dc.size();
}
const V dqs = subset(dx, Span(np*nw, 1, varstart));
varstart += dqs.size();
@@ -1533,8 +1532,6 @@ namespace {
//Polymer concentration updates.
if (has_polymer_) {
const V dc = subset(dx, Span(nc));
int varstart = nc;
const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc, 1);
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, state.concentration().begin());

View File

@@ -85,7 +85,7 @@ namespace Opm
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr<PolymerInflowInterface> polymer_inflow,
const PolymerInflowInterface& polymer_inflow,
const double* gravity,
const bool disgas,
const bool vapoil,

View File

@@ -82,7 +82,7 @@ namespace Opm
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr<PolymerInflowInterface> polymer_inflow,
const PolymerInflowInterface& polymer_inflow,
const double* gravity,
bool has_disgas,
bool has_vapoil,
@@ -116,7 +116,7 @@ namespace Opm
// Solvers
const DerivedGeology& geo_;
NewtonIterationBlackoilInterface& solver_;
std::shared_ptr<PolymerInflowInterface> polymer_inflow_;
const PolymerInflowInterface& polymer_inflow_;
// Misc. data
std::vector<int> allcells_;
const bool has_disgas_;
@@ -148,7 +148,7 @@ namespace Opm
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr<PolymerInflowInterface> polymer_inflow,
const PolymerInflowInterface& polymer_inflow,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
@@ -241,7 +241,7 @@ namespace Opm
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr<PolymerInflowInterface> polymer_inflow,
const PolymerInflowInterface& polymer_inflow,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
@@ -327,7 +327,7 @@ namespace Opm
props_.permeability());
const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoil well_state;
well_state.init(wells, state);
well_state.init(wells, state.blackoilState());
if (timer.currentStepNum() != 0) {
// Transfer previous well state to current.
well_state.partialCopy(prev_well_state, *wells, prev_well_state.numWells());
@@ -345,7 +345,7 @@ namespace Opm
if (timer.currentStepNum() == 0) {
output_writer_.writeInit(timer);
}
output_writer_.writeTimeStep(timer, state, well_state.basicWellState());
output_writer_.writeTimeStep(timer, state.blackoilState(), well_state.basicWellState());
}
// Max oil saturation (for VPPARS), hysteresis update.
@@ -358,10 +358,9 @@ namespace Opm
// compute polymer inflow
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_));
if (has_polymer_) {
polymer_inflow_.reset(new PolymerInflowFromDeck(deck, wells, Opm::UgGridHelpers::numCells(grid_)));
*polymer_inflow_.getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
polymer_inflow_.getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
}
// Run a single step of the solver.
solver_timer.start();
@@ -395,7 +394,7 @@ namespace Opm
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_);
output_writer_.writeTimeStep(timer, state, prev_well_state.basicWellState());
output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState());
}
// Stop timer and create timing report
@@ -539,7 +538,7 @@ namespace Opm
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x);
rateConverter_.defineState(x.blackoilState());
std::vector<double> distr (np);
std::vector<double> hrates(np);