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

Merge pull request #487 from totto82/solventModel

Browse Source 
Changes to the solvent model
This commit is contained in:
Atgeirr Flø Rasmussen 2015-10-06 11:19:58 +02:00
commit 68caa026d9
8 changed files with 305 additions and 338 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -777,8 +777,15 @@ namespace detail {
const V depth = cellCentroidsZToEigen(grid_); const V depth = cellCentroidsZToEigen(grid_);
const V pdepth = subset(depth, well_cells); const V pdepth = subset(depth, well_cells);
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf); std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
// Surface density. // Surface density.
std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases); DataBlock surf_dens(nperf, pu.num_phases);
for (int phase = 0; phase < pu.num_phases; ++ phase) {
surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]);
}
std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
// Gravity // Gravity
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
@ -786,7 +793,7 @@ namespace detail {
std::vector<double> cd = std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities( WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(), wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens); b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 3. Compute pressure deltas // 3. Compute pressure deltas
std::vector<double> cdp = std::vector<double> cdp =
@ -941,7 +948,7 @@ namespace detail {
if (param_.update_equations_scaling_) { if (param_.update_equations_scaling_) {
updateEquationsScaling(); asImpl().updateEquationsScaling();
} }
} }
@ -1538,7 +1545,7 @@ namespace detail {
std::vector<ADB::M> old_derivs = state.qs.derivative(); std::vector<ADB::M> old_derivs = state.qs.derivative();
state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); state.qs = ADB::function(std::move(new_qs), std::move(old_derivs));
} }
computeWellConnectionPressures(state, well_state); asImpl().computeWellConnectionPressures(state, well_state);
} }
} }

40
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -85,12 +85,6 @@ namespace Opm {
ReservoirState& reservoir_state, ReservoirState& reservoir_state,
WellState& well_state); WellState& well_state);
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
/// \param[in] dt timestep length
/// \param[in] iteration current iteration number
bool getConvergence(const double dt, const int iteration);
/// Assemble the residual and Jacobian of the nonlinear system. /// Assemble the residual and Jacobian of the nonlinear system.
/// \param[in] reservoir_state reservoir state variables /// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables /// \param[in, out] well_state well state variables
@ -195,6 +189,11 @@ namespace Opm {
const SolutionState& state, const SolutionState& state,
WellState& xw); WellState& xw);
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
void updateEquationsScaling();
void void
computeMassFlux(const int actph , computeMassFlux(const int actph ,
const V& transi, const V& transi,
@ -205,35 +204,6 @@ namespace Opm {
const std::vector<PhasePresence> const std::vector<PhasePresence>
phaseCondition() const {return this->phaseCondition_;} phaseCondition() const {return this->phaseCondition_;}
/// \brief Compute the reduction within the convergence check.
/// \param[in] B A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. B.col(i) contains the values
/// for phase i.
/// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. tempV.col(i) contains the
/// values
/// for phase i.
/// \param[in] R A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. B.col(i) contains the values
/// for phase i.
/// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
/// of R for the phase i.
/// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
/// maximum of tempV for the phase i.
/// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
/// of B for the phase i.
/// \param[in] nc The number of cells of the local grid.
/// \return The total pore volume over all cells.
double
convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
std::array<double,MaxNumPhases+1>& R_sum,
std::array<double,MaxNumPhases+1>& maxCoeff,
std::array<double,MaxNumPhases+1>& B_avg,
std::vector<double>& maxNormWell,
int nc,
int nw) const;
}; };

410
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -94,10 +94,14 @@ namespace Opm {
// If deck has solvent, residual_ should contain solvent equation. // If deck has solvent, residual_ should contain solvent equation.
rq_.resize(fluid_.numPhases() + 1); rq_.resize(fluid_.numPhases() + 1);
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null()); residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
Base::material_name_.push_back("Solvent");
assert(solvent_pos_ == fluid_.numPhases()); assert(solvent_pos_ == fluid_.numPhases());
if (has_vapoil_) { if (has_vapoil_) {
OPM_THROW(std::runtime_error, "Solvent option only works with dead gas\n"); OPM_THROW(std::runtime_error, "Solvent option only works with dead gas\n");
} }
residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
} }
} }
@ -224,9 +228,18 @@ namespace Opm {
} }
template <class Grid>
void
BlackoilSolventModel<Grid>::updateEquationsScaling()
{
Base::updateEquationsScaling();
assert(MaxNumPhases + 1 == residual_.matbalscale.size());
if (has_solvent_) {
const ADB& temp_b = rq_[solvent_pos_].b;
ADB::V B = 1. / temp_b.value();
residual_.matbalscale[solvent_pos_] = B.mean();
}
}
template <class Grid> template <class Grid>
void BlackoilSolventModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s, void BlackoilSolventModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
@ -250,17 +263,19 @@ namespace Opm {
? state.saturation[ pu.phase_pos[ Gas ] ] ? state.saturation[ pu.phase_pos[ Gas ] ]
: zero); : zero);
Selector<double> zero_selector(ss.value(), Selector<double>::Zero);
V F_solvent = zero_selector.select(ss, ss / (ss + sg)).value();
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf); const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
const int nw = wells().number_of_wells; Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
V wellSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf); ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
const int nw = wells().number_of_wells;
V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
V isProducer = V::Zero(nperf);
V ones = V::Constant(nperf,1.0);
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
if(wells().type[w] == PRODUCER) { if(wells().type[w] == PRODUCER) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
wellSolventFraction[perf] = F_solvent[well_cells[perf]]; isProducer[perf] = 1;
} }
} }
} }
@ -271,7 +286,7 @@ namespace Opm {
// remove contribution from the dissolved gas. // remove contribution from the dissolved gas.
// TODO compensate for gas in the oil phase // TODO compensate for gas in the oil phase
assert(!has_vapoil_); assert(!has_vapoil_);
const ADB cq_s_solvent = wellSolventFraction * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]); const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]);
// Solvent contribution to the mass balance equation is given as a fraction // Solvent contribution to the mass balance equation is given as a fraction
// of the gas contribution. // of the gas contribution.
@ -284,6 +299,142 @@ namespace Opm {
} }
} }
template <class Grid>
void BlackoilSolventModel<Grid>::computeWellConnectionPressures(const SolutionState& state,
const WellState& xw)
{
if( ! Base::localWellsActive() ) return ;
using namespace Opm::AutoDiffGrid;
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Surface density.
const PhaseUsage& pu = fluid_.phaseUsage();
//std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
DataBlock surf_dens(nperf, pu.num_phases);
for (int phase = 0; phase < pu.num_phases; ++ phase) {
surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]);
}
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
DataBlock b(nperf, pu.num_phases);
std::vector<double> rsmax_perf(nperf, 0.0);
std::vector<double> rvmax_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert(active_[Oil]);
const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
if (has_solvent_) {
const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
// A weighted sum of the b-factors of gas and solvent are used.
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB zero = ADB::constant(V::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
V F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
const int nw = wells().number_of_wells;
V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
V isProducer = V::Zero(nperf);
V ones = V::Constant(nperf,1.0);
for (int w = 0; w < nw; ++w) {
if(wells().type[w] == PRODUCER) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
isProducer[perf] = 1;
}
}
}
F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
bg = bg * (ones - F_solvent);
bg = bg + F_solvent * bs;
const V& rhog = surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]);
const V& rhos = solvent_props_.solventSurfaceDensity(well_cells);
surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]) = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
}
// b and surf_dens_perf is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
// Extract well connection depths.
const V depth = cellCentroidsZToEigen(grid_);
const V pdepth = subset(depth, well_cells);
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
// Gravity
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 3. Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_depth, cd, grav);
// 4. Store the results
Base::well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
Base::well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
}
@ -364,7 +515,7 @@ namespace Opm {
? state.saturation[ pu.phase_pos[ Gas ] ] ? state.saturation[ pu.phase_pos[ Gas ] ]
: zero); : zero);
Selector<double> zero_selector(ss.value(), Selector<double>::Zero); Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
ADB F_solvent = zero_selector.select(ss, ss / (ss + sg)); ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
V ones = V::Constant(nc, 1.0); V ones = V::Constant(nc, 1.0);
@ -372,16 +523,23 @@ namespace Opm {
const ADB mu = solvent_props_.muSolvent(phasePressure,cells_); const ADB mu = solvent_props_.muSolvent(phasePressure,cells_);
rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu; rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu;
rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob;
const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b; const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b;
const ADB rhoavg_solvent = ops_.caver * rho_solvent; const ADB rhoavg_solvent = ops_.caver * rho_solvent;
rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix())); rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix()));
UpwindSelector<double> upwind(grid_, ops_, rq_[solvent_pos_].dh.value()); UpwindSelector<double> upwind_solvent(grid_, ops_, rq_[solvent_pos_].dh.value());
// Compute solvent flux. // Compute solvent flux.
rq_[solvent_pos_].mflux = upwind.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh); rq_[solvent_pos_].mflux = upwind_solvent.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh);
// Update gas mobility and flux
rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob;
const ADB& b = rq_[ actph ].b;
const ADB& mob = rq_[ actph ].mob;
const ADB& dh = rq_[ actph ].dh;
UpwindSelector<double> upwind_gas(grid_, ops_, dh.value());
rq_[ actph ].mflux = upwind_gas.select(b * mob) * (transi * dh);
} }
} }
@ -474,10 +632,27 @@ namespace Opm {
// Gas and solvent is combinded and solved together // Gas and solvent is combinded and solved together
// The input in the well equation is then the // The input in the well equation is then the
// total gas phase = hydro carbon gas + solvent gas // total gas phase = hydro carbon gas + solvent gas
// This may need to be reconsidered later, as the model
// is tested. // The total mobility is the sum of the solvent and gas mobiliy
mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells); mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells);
b_perfcells[gas_pos] += subset(rq_[solvent_pos_].b, well_cells);
// A weighted sum of the b-factors of gas and solvent are used.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB zero = ADB::constant(V::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
V ones = V::Constant(nperf,1.0);
b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells));
} }
if (param_.solve_welleq_initially_ && initial_assembly) { if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step // solve the well equations as a pre-processing step
@ -493,207 +668,6 @@ namespace Opm {
} }
template <class Grid>
double
BlackoilSolventModel<Grid>::convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
std::array<double,MaxNumPhases+1>& R_sum,
std::array<double,MaxNumPhases+1>& maxCoeff,
std::array<double,MaxNumPhases+1>& B_avg,
std::vector<double>& maxNormWell,
int nc,
int nw) const
{
// Do the global reductions
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
// Compute the global number of cells and porevolume
std::vector<int> v(nc, 1);
auto nc_and_pv = std::tuple<int, double>(0, 0.0);
auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
for ( int idx=0; idx<MaxNumPhases+1; ++idx )
{
if ((idx == MaxNumPhases && has_solvent_) || active_[idx]) { // Dealing with solvent *or* an active phase.
auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
auto containers = std::make_tuple(B.col(idx),
tempV.col(idx),
R.col(idx));
auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
Opm::Reduction::makeGlobalMaxFunctor<double>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
info.computeReduction(containers, operators, values);
B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
maxCoeff[idx] = std::get<1>(values);
R_sum[idx] = std::get<2>(values);
if (idx != MaxNumPhases) { // We do not compute a well flux residual for solvent.
maxNormWell[idx] = 0.0;
for ( int w=0; w<nw; ++w ) {
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
}
}
}
else
{
maxNormWell[idx] = R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
}
}
info.communicator().max(&maxNormWell[0], MaxNumPhases+1);
// Compute pore volume
return std::get<1>(nc_and_pv);
}
else
#endif
{
for ( int idx=0; idx<MaxNumPhases+1; ++idx )
{
if (((idx == MaxNumPhases && has_solvent_) || active_[idx]) ) { // Dealing with solvent *or* an active phase.
B_avg[idx] = B.col(idx).sum()/nc;
maxCoeff[idx] = tempV.col(idx).maxCoeff();
R_sum[idx] = R.col(idx).sum();
}
else
{
R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
}
if (idx != MaxNumPhases) { // We do not compute a well flux residual for polymer.
maxNormWell[idx] = 0.0;
for ( int w=0; w<nw; ++w ) {
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
}
}
}
// Compute total pore volume
return geo_.poreVolume().sum();
}
}
template <class Grid>
bool
BlackoilSolventModel<Grid>::getConvergence(const double dt, const int iteration)
{
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int nw = wellsActive() ? wells().number_of_wells : 0;
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const V pv = geo_.poreVolume();
const std::vector<PhasePresence> cond = phaseCondition();
std::array<double,MaxNumPhases+1> CNV = {{0., 0., 0., 0.}};
std::array<double,MaxNumPhases+1> R_sum = {{0., 0., 0., 0.}};
std::array<double,MaxNumPhases+1> B_avg = {{0., 0., 0., 0.}};
std::array<double,MaxNumPhases+1> maxCoeff = {{0., 0., 0., 0.}};
std::array<double,MaxNumPhases+1> mass_balance_residual = {{0., 0., 0., 0.}};
std::array<double,MaxNumPhases> well_flux_residual = {{0., 0., 0.}};
std::size_t cols = MaxNumPhases+1; // needed to pass the correct type to Eigen
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> B(nc, cols);
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> R(nc, cols);
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> tempV(nc, cols);
std::vector<double> maxNormWell(MaxNumPhases);
for ( int idx=0; idx<MaxNumPhases; ++idx )
{
if (active_[idx]) {
const int pos = pu.phase_pos[idx];
const ADB& tempB = rq_[pos].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;
}
}
if (has_solvent_) {
const ADB& tempB = rq_[solvent_pos_].b;
B.col(MaxNumPhases) = 1. / tempB.value();
R.col(MaxNumPhases) = residual_.material_balance_eq[solvent_pos_].value();
tempV.col(MaxNumPhases) = R.col(MaxNumPhases).abs()/pv;
}
const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg,
maxNormWell, nc, nw);
bool converged_MB = true;
bool converged_CNV = true;
bool converged_Well = true;
// Finish computation
for ( int idx = 0; idx < (MaxNumPhases + 1) ; ++idx )
{
CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
if (idx != MaxNumPhases) { // No well flux residual for polymer.
well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
}
}
const double residualWell = detail::infinityNormWell(residual_.well_eq,
linsolver_.parallelInformation());
converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
const bool converged = converged_MB && converged_CNV && converged_Well;
// if one of the residuals is NaN, throw exception, so that the solver can be restarted
if (std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Oil]) || mass_balance_residual[Oil] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Gas]) || mass_balance_residual[Gas] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Gas]) || mass_balance_residual[MaxNumPhases] > maxResidualAllowed() ||
std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() ||
std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() ||
std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() ||
std::isnan(CNV[MaxNumPhases]) || CNV[MaxNumPhases] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ||
std::isnan(residualWell) || residualWell > maxResidualAllowed() )
{
OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!");
}
if ( terminal_output_ )
{
// Only rank 0 does print to std::cout
if (iteration == 0) {
std::cout << "\nIter MB(WATER) MB(OIL) MB(GAS) MB(SOLVENT) CNVW CNVO CNVG CNVS W-FLUX(W) W-FLUX(O) W-FLUX(G)\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[MaxNumPhases]
<< std::setw(11) << CNV[Water]
<< std::setw(11) << CNV[Oil]
<< std::setw(11) << CNV[Gas]
<< std::setw(11) << CNV[MaxNumPhases]
<< std::setw(11) << well_flux_residual[Water]
<< std::setw(11) << well_flux_residual[Oil]
<< std::setw(11) << well_flux_residual[Gas]
<< std::endl;
std::cout.precision(oprec);
std::cout.flags(oflags);
}
return converged;
}
} }

8
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
View File

@ -97,9 +97,10 @@ namespace Opm
const Wells* wells) const Wells* wells)
{ {
// compute solvent inflow // compute solvent inflow
const int nw = wells->number_of_wells;
std::vector<double> perfcells_fraction(wells->well_connpos[nw], 0.0);
if (deck_->hasKeyword("WSOLVENT")) { if (deck_->hasKeyword("WSOLVENT")) {
const int nw = wells->number_of_wells;
std::vector<double> perfcells_fraction(wells->well_connpos[nw]);
size_t currentStep = timer.currentStepNum(); size_t currentStep = timer.currentStepNum();
ScheduleConstPtr schedule = BaseType::eclipse_state_->getSchedule(); ScheduleConstPtr schedule = BaseType::eclipse_state_->getSchedule();
@ -135,11 +136,10 @@ namespace Opm
} }
} }
} }
well_state.solventFraction() = perfcells_fraction;
} }
well_state.solventFraction() = perfcells_fraction;
} }
} // namespace Opm } // namespace Opm
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_IMPL_HEADER_INCLUDED #endif // OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_IMPL_HEADER_INCLUDED

144
opm/autodiff/SolventPropsAdFromDeck.cpp
View File

@ -36,87 +36,89 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
const int number_of_cells, const int number_of_cells,
const int* global_cell) const int* global_cell)
{ {
// retrieve the cell specific PVT table index from the deck if (deck->hasKeyword("SOLVENT")) {
// and using the grid... // retrieve the cell specific PVT table index from the deck
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell); // and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, deck, number_of_cells, global_cell);
// surface densities // surface densities
if (deck->hasKeyword("SDENSITY")) { if (deck->hasKeyword("SDENSITY")) {
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("SDENSITY"); Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("SDENSITY");
int numRegions = densityKeyword->size(); int numRegions = densityKeyword->size();
solvent_surface_densities_.resize(numRegions); solvent_surface_densities_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
solvent_surface_densities_[regionIdx] solvent_surface_densities_[regionIdx]
= densityKeyword->getRecord(regionIdx)->getItem("SOLVENT_DENSITY")->getSIDouble(0); = densityKeyword->getRecord(regionIdx)->getItem("SOLVENT_DENSITY")->getSIDouble(0);
}
} else {
OPM_THROW(std::runtime_error, "SDENSITY must be specified in SOLVENT runs\n");
} }
} else {
OPM_THROW(std::runtime_error, "SDENSITY must be specified in SOLVENT runs\n");
}
auto tables = eclState->getTableManager(); auto tables = eclState->getTableManager();
// pvt // pvt
const TableContainer& pvdsTables = tables->getPvdsTables(); const TableContainer& pvdsTables = tables->getPvdsTables();
if (!pvdsTables.empty()) { if (!pvdsTables.empty()) {
int numRegions = pvdsTables.size(); int numRegions = pvdsTables.size();
// resize the attributes of the object // resize the attributes of the object
b_.resize(numRegions); b_.resize(numRegions);
viscosity_.resize(numRegions); viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions); inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdsTable& pvdsTable = pvdsTables.getTable<PvdsTable>(regionIdx); const Opm::PvdsTable& pvdsTable = pvdsTables.getTable<PvdsTable>(regionIdx);
// Copy data // Copy data
const std::vector<double>& press = pvdsTable.getPressureColumn(); const std::vector<double>& press = pvdsTable.getPressureColumn();
const std::vector<double>& b = pvdsTable.getFormationFactorColumn(); const std::vector<double>& b = pvdsTable.getFormationFactorColumn();
const std::vector<double>& visc = pvdsTable.getViscosityColumn(); const std::vector<double>& visc = pvdsTable.getViscosityColumn();
const int sz = b.size(); const int sz = b.size();
std::vector<double> inverseB(sz); std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) { for (int i = 0; i < sz; ++i) {
inverseB[i] = 1.0 / b[i]; inverseB[i] = 1.0 / b[i];
}
std::vector<double> inverseBmu(sz);
for (int i = 0; i < sz; ++i) {
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
} else {
OPM_THROW(std::runtime_error, "PVDS must be specified in SOLVENT runs\n");
}
const TableContainer& ssfnTables = tables->getSsfnTables();
// relative permeabilty multiplier
if (!ssfnTables.empty()) {
int numRegions = pvdsTables.size();
if(numRegions > 1) {
OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN");
}
// resize the attributes of the object
krg_.resize(numRegions);
krs_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::SsfnTable& ssfnTable = ssfnTables.getTable<SsfnTable>(regionIdx);
// Copy data
const std::vector<double>& solventFraction = ssfnTable.getSolventFractionColumn();
const std::vector<double>& krg = ssfnTable.getGasRelPermMultiplierColumn();
const std::vector<double>& krs = ssfnTable.getSolventRelPermMultiplierColumn();
krg_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krg);
krs_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krs);
} }
std::vector<double> inverseBmu(sz); } else {
for (int i = 0; i < sz; ++i) { OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n");
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
} }
} else {
OPM_THROW(std::runtime_error, "PVDS must be specified in SOLVENT runs\n");
}
const TableContainer& ssfnTables = tables->getSsfnTables();
// relative permeabilty multiplier
if (!ssfnTables.empty()) {
int numRegions = pvdsTables.size();
if(numRegions > 1) {
OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN");
}
// resize the attributes of the object
krg_.resize(numRegions);
krs_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::SsfnTable& ssfnTable = ssfnTables.getTable<SsfnTable>(regionIdx);
// Copy data
const std::vector<double>& solventFraction = ssfnTable.getSolventFractionColumn();
const std::vector<double>& krg = ssfnTable.getGasRelPermMultiplierColumn();
const std::vector<double>& krs = ssfnTable.getSolventRelPermMultiplierColumn();
krg_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krg);
krs_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krs);
}
} else {
OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n");
} }
} }

11
opm/autodiff/WellDensitySegmented.cpp
View File

@ -36,7 +36,7 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf, const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens) const std::vector<double>& surf_dens_perf)
{ {
// Verify that we have consistent input. // Verify that we have consistent input.
const int np = wells.number_of_phases; const int np = wells.number_of_phases;
@ -45,8 +45,8 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
if (wells.number_of_phases != phase_usage.num_phases) { if (wells.number_of_phases != phase_usage.num_phases) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. phase_usage."); OPM_THROW(std::logic_error, "Inconsistent input: wells vs. phase_usage.");
} }
if (surf_dens.size() != size_t(wells.number_of_phases)) { if (nperf*np != int(surf_dens_perf.size())) {
OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage."); OPM_THROW(std::logic_error, "Inconsistent input: wells vs. surf_dens.");
} }
if (nperf*np != int(wstate.perfPhaseRates().size())) { if (nperf*np != int(wstate.perfPhaseRates().size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate."); OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
@ -94,6 +94,7 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid]; const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
std::vector<double> mix(np); std::vector<double> mix(np);
std::vector<double> x(np); std::vector<double> x(np);
std::vector<double> surf_dens(np);
std::vector<double> dens(nperf); std::vector<double> dens(nperf);
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) { for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
@ -130,6 +131,10 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells,
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
volrat += x[phase] / b_perf[perf*np + phase]; volrat += x[phase] / b_perf[perf*np + phase];
} }
for (int phase = 0; phase < np; ++phase) {
surf_dens[phase] = surf_dens_perf[perf*np + phase];
}
// Compute segment density. // Compute segment density.
dens[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; dens[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
} }

4
opm/autodiff/WellDensitySegmented.hpp
View File

@ -47,14 +47,14 @@ namespace Opm
/// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation /// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation
/// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N /// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N
/// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N /// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N
/// \param[in] surf_dens surface densities for active components, size P /// \param[in] surf_dens surface densities for active components, size NP, P values per perforation
static std::vector<double> computeConnectionDensities(const Wells& wells, static std::vector<double> computeConnectionDensities(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate, const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage, const PhaseUsage& phase_usage,
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf, const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens); const std::vector<double>& surf_dens_perf);
/// Compute pressure deltas. /// Compute pressure deltas.
/// Notation: N = number of perforations, P = number of phases. /// Notation: N = number of perforations, P = number of phases.

11
tests/test_welldensitysegmented.cpp
View File

@ -88,7 +88,16 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
const std::vector<double> rsmax_perf = { 50, 50, 50, 50, 50, 50, 50, 50, 50, 50 }; const std::vector<double> rsmax_perf = { 50, 50, 50, 50, 50, 50, 50, 50, 50, 50 };
const std::vector<double> rvmax_perf = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 }; const std::vector<double> rvmax_perf = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
const std::vector<double> z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 }; const std::vector<double> z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 };
const std::vector<double> surf_dens = { 1000.0, 800.0, 10.0 }; const std::vector<double> surf_dens = { 1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0,
1000.0, 800.0, 10.0};
const double gravity = Opm::unit::gravity; const double gravity = Opm::unit::gravity;
std::vector<double> cd = std::vector<double> cd =