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

Enable single-phase runs.

Browse Source
This commit is contained in:
Halvor M. Nilsen 2019-10-09 15:24:23 +02:00 committed by Atgeirr Flø Rasmussen
parent e4e8425bad
commit d873ae165d
12 changed files with 264 additions and 171 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
ebos/ecloutputblackoilmodule.hh
View File

@ -395,7 +395,16 @@ public:
}
if (oilPressure_.size() > 0) {
oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx));
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx));
}else{
// put pressure in oil pressure for output
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(waterPhaseIdx));
} else {
oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(gasPhaseIdx));
}
}
Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]);
}

17
ebos/eclproblem.hh
View File

@ -665,17 +665,20 @@ public:
minTimeStepSize_ = tuning.getTSMINZ(0);
}
initFluidSystem_();
// deal with DRSDT
unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
maxDRs_.resize(ntpvt, 1e30);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
size_t numDof = this->model().numGridDof();
lastRs_.resize(numDof, 0.0);
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxOilSaturation_.resize(numDof, 0.0);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
maxDRs_.resize(ntpvt, 1e30);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
lastRs_.resize(numDof, 0.0);
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxOilSaturation_.resize(numDof, 0.0);
}
initFluidSystem_();
updateElementDepths_();
readRockParameters_();
readMaterialParameters_();

7
opm/core/props/phaseUsageFromDeck.hpp
View File

@ -54,14 +54,9 @@ namespace Opm
}
}
// Only 2 or 3 phase systems handled.
if (pu.num_phases < 2 || pu.num_phases > 3) {
OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
}
// We need oil systems, since we do not support the keywords needed for
// water-gas systems.
if (!pu.phase_used[BlackoilPhases::Liquid]) {
if (!pu.phase_used[BlackoilPhases::Liquid] && !(pu.num_phases == 1)) {
OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
}

43
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@ -422,29 +422,34 @@ namespace Opm {
Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
Scalar oilSaturationOld = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
}
if (FluidSystem::numActivePhases() > 1) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
{
saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
}
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
resultDelta += tmp*tmp;
resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
assert(std::isfinite(resultDelta));
assert(std::isfinite(resultDenom));
}
}
// NB fix me! adding pressures changes to satutation changes does not make sense
Scalar tmp = pressureNew - pressureOld;
resultDelta += tmp*tmp;
resultDenom += pressureNew*pressureNew;
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
const Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
resultDelta += tmpSat * tmpSat;
resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
}
}
resultDelta = gridView.comm().sum(resultDelta);

6
opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
View File

@ -60,6 +60,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
{
++current_step_;
current_time_ += dt_;
assert(dt_ > 0);
// store used time step sizes
steps_.push_back( dt_ );
return *this;
@ -71,7 +72,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
double remaining = (total_time_ - current_time_);
// apply max time step if it was set
dt_ = std::min( dt_estimate, max_time_step_ );
assert(dt_ > 0);
if( remaining > 0 ) {
// set new time step (depending on remaining time)
@ -81,6 +82,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
if( dt_ > max_time_step_ ) {
dt_ = 0.5 * remaining;
}
assert(dt_ > 0);
return;
}
@ -89,6 +91,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
if( 1.5 * dt_ > remaining ) {
dt_ = 0.5 * remaining;
assert(dt_ > 0);
return;
}
}
@ -102,6 +105,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
double AdaptiveSimulatorTimer::currentStepLength () const
{
assert(dt_ > 0);
return dt_;
}

5
opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp
View File

@ -293,9 +293,10 @@ namespace Opm {
double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
substepTimer.simulationTimeElapsed());
assert(dtEstimate > 0);
// limit the growth of the timestep size by the growth factor
dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
assert(dtEstimate > 0);
// further restrict time step size growth after convergence problems
if (restarts > 0) {
dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
@ -309,7 +310,7 @@ namespace Opm {
const double maxPredictionTHPTimestep = 16.0 * unit::day;
dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
}
assert(dtEstimate > 0);
if (timestepVerbose_) {
std::ostringstream ss;
substepReport.reportStep(ss);

8
opm/simulators/timestepping/TimeStepControl.cpp
View File

@ -128,13 +128,17 @@ namespace Opm
{
// shift errors
for( int i=0; i<2; ++i ) {
errors_[ i ] = errors_[i+1];
errors_[ i ] = errors_[i+1];
}
// store new error
const double error = relChange.relativeChange();
errors_[ 2 ] = error;
for( int i=0; i<2; ++i ) {
assert(std::isfinite(errors_[i]));
assert(errors_[i]>0);
}
if( error > tol_ )
{
// adjust dt by given tolerance

20
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -251,9 +251,19 @@ namespace Opm {
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
cellPressures[cellIdx] = p;
// copy of get perfpressure in Standard well
// exept for value
double perf_pressure = 0.0;
if (Indices::oilEnabled) {
perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
} else {
if (Indices::waterEnabled) {
perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
} else {
perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
}
}
cellPressures[cellIdx] = perf_pressure;
}
well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_);
@ -1525,8 +1535,8 @@ namespace Opm {
int
BlackoilWellModel<TypeTag>::numComponents() const
{
if (numPhases() == 2) {
return 2;
if (numWells() > 0 && numPhases() < 3) {
return numPhases();
}
int numComp = FluidSystem::numComponents;
if (has_solvent_) {

12
opm/simulators/wells/RateConverter.hpp
View File

@ -493,11 +493,13 @@ namespace Opm {
// sum p, rs, rv, and T.
double hydrocarbonPV = pv_cell*hydrocarbon;
pv += hydrocarbonPV;
p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
rs += fs.Rs().value()*hydrocarbonPV;
rv += fs.Rv().value()*hydrocarbonPV;
T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
if (hydrocarbonPV > 0) {
pv += hydrocarbonPV;
p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
rs += fs.Rs().value()*hydrocarbonPV;
rv += fs.Rv().value()*hydrocarbonPV;
T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
}
}
for (const auto& reg : rmap_.activeRegions()) {

3
opm/simulators/wells/StandardWell.hpp
View File

@ -59,6 +59,7 @@ namespace Opm
using typename Base::Indices;
using typename Base::RateConverterType;
using typename Base::SparseMatrixAdapter;
using typename Base::FluidState;
using Base::numEq;
@ -293,6 +294,8 @@ namespace Opm
EvalWell extendEval(const Eval& in) const;
Eval getPerfPressure(const FluidState& fs) const;
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;

294
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -238,6 +238,10 @@ namespace Opm
StandardWell<TypeTag>::
wellVolumeFraction(const unsigned compIdx) const
{
if (FluidSystem::numActivePhases() == 1) {
return EvalWell(numWellEq_ + numEq, 1.0);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[WFrac];
}
@ -305,6 +309,24 @@ namespace Opm
template<typename TypeTag>
typename StandardWell<TypeTag>::Eval
StandardWell<TypeTag>::getPerfPressure(const typename StandardWell<TypeTag>::FluidState& fs) const
{
Eval pressure;
if (Indices::oilEnabled) {
pressure = fs.pressure(FluidSystem::oilPhaseIdx);
} else {
if (Indices::waterEnabled) {
pressure = fs.pressure(FluidSystem::waterPhaseIdx);
} else {
pressure = fs.pressure(FluidSystem::gasPhaseIdx);
}
}
return pressure;
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
@ -321,7 +343,7 @@ namespace Opm
{
const auto& fs = intQuants.fluidState();
const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const EvalWell pressure = extendEval(getPerfPressure(fs));
const EvalWell rs = extendEval(fs.Rs());
const EvalWell rv = extendEval(fs.Rv());
std::vector<EvalWell> b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0});
@ -669,7 +691,9 @@ namespace Opm
|| (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) {
const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value();
const auto& fs = intQuants.fluidState();
Eval perf_pressure = getPerfPressure(fs);
const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value();
const bool new_well = schedule.hasWellEvent(name(), ScheduleEvents::NEW_WELL, current_step_);
double productivity_index = cq_s[compIdx].value() / drawdown;
scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger);
@ -679,12 +703,15 @@ namespace Opm
}
// add vol * dF/dt + Q to the well equations;
for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
// TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
// since all the rates are under surface condition
EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
EvalWell resWell_loc(numWellEq_ + numEq, 0.0);
if (FluidSystem::numActivePhases() > 1) {
assert(dt > 0);
resWell_loc += (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
}
resWell_loc -= getQs(componentIdx) * well_efficiency_factor_;
for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
@ -919,29 +946,30 @@ namespace Opm
const double relaxation_factor_fractions = (well_type_ == PRODUCER) ?
relaxationFactorFractionsProducer(old_primary_variables, dwells)
: 1.0;
if (FluidSystem::numActivePhases() > 1) {
// update the second and third well variable (The flux fractions)
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
// primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
}
// update the second and third well variable (The flux fractions)
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
// primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
}
if (has_solvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
}
processFractions();
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
}
if (has_solvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
}
processFractions();
// updating the total rates Q_t
const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells);
primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
@ -956,6 +984,10 @@ namespace Opm
}
updateExtraPrimaryVariables(dwells);
for (int i = 0; i < primary_variables_.size(); ++i) {
assert(Opm::isfinite(primary_variables_[i]));
}
}
@ -1073,49 +1105,53 @@ namespace Opm
updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const
{
const PhaseUsage& pu = phaseUsage();
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil];
std::vector<double> F(number_of_phases_, 0.0);
F[oil_pos] = 1.0;
if (FluidSystem::numActivePhases() == 1) {
F[0] = 1.0;
} else {
assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
const int oil_pos = pu.phase_pos[Oil];
F[oil_pos] = 1.0;
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = primary_variables_[WFrac];
F[oil_pos] -= F[water_pos];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = primary_variables_[GFrac];
F[oil_pos] -= F[gas_pos];
}
double F_solvent = 0.0;
if (has_solvent) {
F_solvent = primary_variables_[SFrac];
F[oil_pos] -= F_solvent;
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < number_of_phases_; ++p) {
const double scal = scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scal > 0) {
F[p] /= scal ;
} else {
// this should only happens to injection wells
F[p] = 0.;
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
const int water_pos = pu.phase_pos[Water];
F[water_pos] = primary_variables_[WFrac];
F[oil_pos] -= F[water_pos];
}
}
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
// More testing is needed to make sure this is correct for well groups and THP.
if (has_solvent){
F_solvent /= scalingFactor(contiSolventEqIdx);
F[pu.phase_pos[Gas]] += F_solvent;
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
const int gas_pos = pu.phase_pos[Gas];
F[gas_pos] = primary_variables_[GFrac];
F[oil_pos] -= F[gas_pos];
}
double F_solvent = 0.0;
if (has_solvent) {
F_solvent = primary_variables_[SFrac];
F[oil_pos] -= F_solvent;
}
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < number_of_phases_; ++p) {
const double scal = scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scal > 0) {
F[p] /= scal ;
} else {
// this should only happens to injection wells
F[p] = 0.;
}
}
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
// More testing is needed to make sure this is correct for well groups and THP.
if (has_solvent){
F_solvent /= scalingFactor(contiSolventEqIdx);
F[pu.phase_pos[Gas]] += F_solvent;
}
}
well_state.bhp()[index_of_well_] = primary_variables_[Bhp];
// calculate the phase rates based on the primary variables
@ -1292,6 +1328,10 @@ namespace Opm
break;
} // end of switch
for (int i = 0; i < primary_variables_.size(); ++i) {
assert(Opm::isfinite(primary_variables_[i]));
}
}
@ -1324,7 +1364,8 @@ namespace Opm
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quantities.fluidState();
// the pressure of the reservoir grid block the well connection is in
const double p_r = fs.pressure(FluidSystem::oilPhaseIdx).value();
Eval perf_pressure = getPerfPressure(fs);
double p_r = perf_pressure.value();
// calculating the b for the connection
std::vector<double> b_perf(num_components_);
@ -1550,8 +1591,9 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
Eval perf_pressure = getPerfPressure(fs);
const double pressure = perf_pressure.value();
const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
const double bhp = getBhp().value();
// Pressure drawdown (also used to determine direction of flow)
@ -1755,6 +1797,7 @@ namespace Opm
// TODO: to check why should be perf - 1
const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1];
const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
// strange choice
const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
if (waterPresent) {
@ -2518,57 +2561,57 @@ namespace Opm
}
}
if (std::abs(total_well_rate) > 0.) {
const auto pu = phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
}
if (has_solvent) {
primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
}
} else { // total_well_rate == 0
if (well_type_ == INJECTOR) {
auto phase = well_ecl_.getInjectionProperties().injectorType;
// only single phase injection handled
if (FluidSystem::numActivePhases() > 1) {
if (std::abs(total_well_rate) > 0.) {
const auto pu = phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (phase == Well2::InjectorType::WATER) {
primary_variables_[WFrac] = 1.0;
} else {
primary_variables_[WFrac] = 0.0;
}
primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (phase == Well2::InjectorType::GAS) {
primary_variables_[GFrac] = 1.0 - wsolvent();
if (has_solvent) {
primary_variables_[SFrac] = wsolvent();
primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
}
if (has_solvent) {
primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
}
} else { // total_well_rate == 0
if (well_type_ == INJECTOR) {
auto phase = well_ecl_.getInjectionProperties().injectorType;
// only single phase injection handled
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
if (phase == Well2::InjectorType::WATER) {
primary_variables_[WFrac] = 1.0;
} else {
primary_variables_[WFrac] = 0.0;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if (phase == Well2::InjectorType::GAS) {
primary_variables_[GFrac] = 1.0 - wsolvent();
if (has_solvent) {
primary_variables_[SFrac] = wsolvent();
}
} else {
primary_variables_[GFrac] = 0.0;
}
} else {
primary_variables_[GFrac] = 0.0;
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type_ == PRODUCER) { // producers
// TODO: the following are not addressed for the solvent case yet
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = 1.0 / np;
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type_ == PRODUCER) { // producers
// TODO: the following are not addressed for the solvent case yet
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = 1.0 / np;
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = 1.0 / np;
}
} else {
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
primary_variables_[GFrac] = 1.0 / np;
}
} else {
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
}
}
// BHP
primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
@ -2579,6 +2622,10 @@ namespace Opm
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
}
}
for (int i = 0; i < primary_variables_.size(); ++i) {
assert(Opm::isfinite(primary_variables_[i]));
}
}
@ -2824,28 +2871,31 @@ namespace Opm
// 0.95 is a experimental value, which remains to be optimized
double relaxation_factor = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
}
if (FluidSystem::numActivePhases() > 1) {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
// not be negative oil fraction later
const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
const double possible_updated_sum = original_sum - relaxed_update;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
&& FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so
// there will not be negative oil fraction later
const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
const double possible_updated_sum = original_sum - relaxed_update;
if (possible_updated_sum > 1.0) {
assert(relaxed_update != 0.);
if (possible_updated_sum > 1.0) {
assert(relaxed_update != 0.);
const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
relaxation_factor *= further_relaxation_factor;
const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
relaxation_factor *= further_relaxation_factor;
}
}
}

9
opm/simulators/wells/WellInterface.hpp
View File

@ -94,6 +94,7 @@ namespace Opm
static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
static const bool has_energy = GET_PROP_VALUE(TypeTag, EnableEnergy);
static const bool has_temperature = GET_PROP_VALUE(TypeTag, EnableTemperature);
// flag for polymer molecular weight related
static const bool has_polymermw = GET_PROP_VALUE(TypeTag, EnablePolymerMW);
static const bool has_foam = GET_PROP_VALUE(TypeTag, EnableFoam);
@ -106,7 +107,13 @@ namespace Opm
// For the conversion between the surface volume rate and reservoir voidage rate
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
static const bool compositionSwitchEnabled = Indices::gasEnabled;
using FluidState = Opm::BlackOilFluidState<Eval,
FluidSystem,
has_temperature,
has_energy,
compositionSwitchEnabled,
Indices::numPhases >;
/// Constructor
WellInterface(const Well2& well, const int time_step, const Wells* wells,
const ModelParameters& param,