mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Fix the 2p simulator
Only tested for oil+water case The blockmatrix and vectors are hardcoded to be 3 and a trivial equation is used for the Gas phase.
This commit is contained in:
Tor Harald Sandve
Andreas Lauser
parent
5e0804b39f
commit
d23270c98f
@ -622,99 +622,102 @@ namespace Opm {
|
||||
double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
|
||||
so -= step * dso;
|
||||
|
||||
// const double drmaxrel = drMaxRel();
|
||||
// Update rs and rv
|
||||
if (has_disgas_) {
|
||||
double& rs = reservoir_state.gasoilratio()[cell_idx];
|
||||
rs -= drs;
|
||||
rs = std::max(rs, 0.0);
|
||||
|
||||
}
|
||||
if (has_vapoil_) {
|
||||
double& rv = reservoir_state.rv()[cell_idx];
|
||||
rv -= drv;
|
||||
rv = std::max(rv, 0.0);
|
||||
}
|
||||
|
||||
// Sg is used as primal variable for water only cells.
|
||||
const double epsilon = 1e-4; //std::sqrt(std::numeric_limits<double>::epsilon());
|
||||
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
|
||||
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
|
||||
double& rs = reservoir_state.gasoilratio()[cell_idx];
|
||||
double& rv = reservoir_state.rv()[cell_idx];
|
||||
|
||||
|
||||
// phase translation sg <-> rs
|
||||
const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
|
||||
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
|
||||
const auto& fs = intQuants.fluidState();
|
||||
switch (hydroCarbonState) {
|
||||
case HydroCarbonState::GasAndOil: {
|
||||
|
||||
if (sw > (1.0 - epsilon)) // water only i.e. do nothing
|
||||
break;
|
||||
|
||||
if (sg <= 0.0 && has_disgas_) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
|
||||
sg = 0;
|
||||
so = 1.0 - sw - sg;
|
||||
const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
// phase for when oil and gas
|
||||
if (active_[Gas] && active_[Oil] ) {
|
||||
// const double drmaxrel = drMaxRel();
|
||||
// Update rs and rv
|
||||
if (has_disgas_) {
|
||||
double& rs = reservoir_state.gasoilratio()[cell_idx];
|
||||
rs = rsSat*(1-epsilon);
|
||||
} else if (so <= 0.0 && has_vapoil_) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv
|
||||
so = 0;
|
||||
sg = 1.0 - sw - so;
|
||||
rs -= drs;
|
||||
rs = std::max(rs, 0.0);
|
||||
|
||||
}
|
||||
if (has_vapoil_) {
|
||||
double& rv = reservoir_state.rv()[cell_idx];
|
||||
// use gas pressure?
|
||||
rv -= drv;
|
||||
rv = std::max(rv, 0.0);
|
||||
}
|
||||
|
||||
// Sg is used as primal variable for water only cells.
|
||||
const double epsilon = 1e-4; //std::sqrt(std::numeric_limits<double>::epsilon());
|
||||
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
|
||||
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
|
||||
double& rs = reservoir_state.gasoilratio()[cell_idx];
|
||||
double& rv = reservoir_state.rv()[cell_idx];
|
||||
|
||||
|
||||
// phase translation sg <-> rs
|
||||
const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
|
||||
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
|
||||
const auto& fs = intQuants.fluidState();
|
||||
switch (hydroCarbonState) {
|
||||
case HydroCarbonState::GasAndOil: {
|
||||
|
||||
if (sw > (1.0 - epsilon)) // water only i.e. do nothing
|
||||
break;
|
||||
|
||||
if (sg <= 0.0 && has_disgas_) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
|
||||
sg = 0;
|
||||
so = 1.0 - sw - sg;
|
||||
const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
double& rs = reservoir_state.gasoilratio()[cell_idx];
|
||||
rs = rsSat*(1-epsilon);
|
||||
} else if (so <= 0.0 && has_vapoil_) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv
|
||||
so = 0;
|
||||
sg = 1.0 - sw - so;
|
||||
double& rv = reservoir_state.rv()[cell_idx];
|
||||
// use gas pressure?
|
||||
const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
rv = rvSat*(1-epsilon);
|
||||
}
|
||||
break;
|
||||
}
|
||||
case HydroCarbonState::OilOnly: {
|
||||
if (sw > (1.0 - epsilon)) {
|
||||
// water only change to Sg
|
||||
rs = 0;
|
||||
rv = 0;
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
|
||||
break;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
if (rs > ( rsSat * (1+epsilon) ) ) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
sg = epsilon;
|
||||
so -= epsilon;
|
||||
rs = rsSat;
|
||||
}
|
||||
break;
|
||||
}
|
||||
case HydroCarbonState::GasOnly: {
|
||||
if (sw > (1.0 - epsilon)) {
|
||||
// water only change to Sg
|
||||
rs = 0;
|
||||
rv = 0;
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
|
||||
break;
|
||||
}
|
||||
const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
rv = rvSat*(1-epsilon);
|
||||
}
|
||||
break;
|
||||
}
|
||||
case HydroCarbonState::OilOnly: {
|
||||
if (sw > (1.0 - epsilon)) {
|
||||
// water only change to Sg
|
||||
rs = 0;
|
||||
rv = 0;
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
|
||||
if (rv > rvSat * (1+epsilon) ) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
so = epsilon;
|
||||
rv = rvSat;
|
||||
sg -= epsilon;
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
if (rs > ( rsSat * (1+epsilon) ) ) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
sg = epsilon;
|
||||
so -= epsilon;
|
||||
rs = rsSat;
|
||||
default:
|
||||
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState);
|
||||
}
|
||||
break;
|
||||
}
|
||||
case HydroCarbonState::GasOnly: {
|
||||
if (sw > (1.0 - epsilon)) {
|
||||
// water only change to Sg
|
||||
rs = 0;
|
||||
rv = 0;
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
|
||||
break;
|
||||
}
|
||||
const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
|
||||
if (rv > rvSat * (1+epsilon) ) {
|
||||
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
|
||||
so = epsilon;
|
||||
rv = rvSat;
|
||||
sg -= epsilon;
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
default:
|
||||
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState);
|
||||
}
|
||||
}
|
||||
|
||||
@ -1237,56 +1240,61 @@ namespace Opm {
|
||||
cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
|
||||
|
||||
// set switching variable and interpretation
|
||||
if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
|
||||
{
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
|
||||
if (active_[Gas] ) {
|
||||
if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
|
||||
{
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
|
||||
}
|
||||
else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
|
||||
{
|
||||
// this case (-> gas only with vaporized oil in the gas) is
|
||||
// relatively expensive as it requires to compute the capillary
|
||||
// pressure in order to get the gas phase pressure. (the reason why
|
||||
// ebos uses the gas pressure here is that it makes the common case
|
||||
// of the primary variable switching code fast because to determine
|
||||
// whether the oil phase appears one needs to compute the Rv value
|
||||
// for the saturated gas phase and if this is not available as a
|
||||
// primary variable, it needs to be computed.) luckily for here, the
|
||||
// gas-only case is not too common, so the performance impact of this
|
||||
// is limited.
|
||||
typedef Opm::SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystem,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false> SatOnlyFluidState;
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
|
||||
|
||||
double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
|
||||
const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
|
||||
double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
|
||||
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
|
||||
}
|
||||
else
|
||||
{
|
||||
assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
|
||||
}
|
||||
} else {
|
||||
// for oil-water case oil pressure should be used as primary variable
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
|
||||
}
|
||||
else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
|
||||
{
|
||||
// this case (-> gas only with vaporized oil in the gas) is
|
||||
// relatively expensive as it requires to compute the capillary
|
||||
// pressure in order to get the gas phase pressure. (the reason why
|
||||
// ebos uses the gas pressure here is that it makes the common case
|
||||
// of the primary variable switching code fast because to determine
|
||||
// whether the oil phase appears one needs to compute the Rv value
|
||||
// for the saturated gas phase and if this is not available as a
|
||||
// primary variable, it needs to be computed.) luckily for here, the
|
||||
// gas-only case is not too common, so the performance impact of this
|
||||
// is limited.
|
||||
typedef Opm::SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystem,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false> SatOnlyFluidState;
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
|
||||
|
||||
double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
|
||||
const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
|
||||
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
|
||||
double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
|
||||
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
|
||||
}
|
||||
else
|
||||
{
|
||||
assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
|
||||
cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
|
||||
cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
|
||||
cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -233,7 +233,7 @@ public:
|
||||
|
||||
const Wells* wells = wells_manager.c_wells();
|
||||
WellState well_state;
|
||||
well_state.init(wells, state, prev_well_state);
|
||||
well_state.init(wells, state, prev_well_state, props_.phaseUsage());
|
||||
|
||||
// give the polymer and surfactant simulators the chance to do their stuff
|
||||
handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
|
||||
|
||||
@ -254,7 +254,7 @@ namespace Opm
|
||||
std::unordered_set<std::string>());
|
||||
|
||||
const Wells* wells = wellsmanager.c_wells();
|
||||
wellstate.resize(wells, simulatorstate); //Resize for restart step
|
||||
wellstate.resize(wells, simulatorstate, phaseusage ); //Resize for restart step
|
||||
auto restarted = Opm::init_from_restart_file(
|
||||
eclipseState_,
|
||||
Opm::UgGridHelpers::numCells(grid) );
|
||||
|
||||
@ -55,6 +55,12 @@
|
||||
|
||||
namespace Opm {
|
||||
|
||||
enum WellVariablePositions {
|
||||
XvarWell = 0,
|
||||
WFrac = 1,
|
||||
GFrac = 2
|
||||
};
|
||||
|
||||
|
||||
/// Class for handling the standard well model.
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
@ -62,15 +68,17 @@ namespace Opm {
|
||||
public:
|
||||
|
||||
// --------- Types ---------
|
||||
typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
|
||||
typedef WellStateFullyImplicitBlackoilDense WellState;
|
||||
typedef BlackoilModelParameters ModelParameters;
|
||||
|
||||
typedef double Scalar;
|
||||
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
|
||||
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
|
||||
static const int blocksize = 3;
|
||||
typedef Dune::FieldVector<Scalar, blocksize > VectorBlockType;
|
||||
typedef Dune::FieldMatrix<Scalar, blocksize, blocksize > MatrixBlockType;
|
||||
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
|
||||
typedef Dune::BlockVector<VectorBlockType> BVector;
|
||||
typedef DenseAd::Evaluation<double, /*size=*/blocksize*2> EvalWell;
|
||||
|
||||
|
||||
// --------- Public methods ---------
|
||||
StandardWellsDense(const Wells* wells_arg,
|
||||
@ -106,6 +114,10 @@ namespace Opm {
|
||||
return;
|
||||
}
|
||||
|
||||
// assumes the gas fractions are stored after water fractions
|
||||
// WellVariablePositions needs to be changed for 2p runs
|
||||
assert (np == 3 || np == 2 && !pu.phase_used[Gas] );
|
||||
|
||||
fluid_ = fluid_arg;
|
||||
active_ = active_arg;
|
||||
vfp_properties_ = vfp_properties_arg;
|
||||
@ -235,10 +247,16 @@ namespace Opm {
|
||||
for (int p2 = 0; p2 < np; ++p2) {
|
||||
if (!only_wells) {
|
||||
ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2);
|
||||
duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].derivative(p2+3); // intput in transformed matrix
|
||||
duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].derivative(p2+blocksize); // intput in transformed matrix
|
||||
duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2);
|
||||
}
|
||||
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+3);
|
||||
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+blocksize);
|
||||
}
|
||||
|
||||
// add trivial equation for 2p cases (Only support water + oil)
|
||||
if (np == 2) {
|
||||
assert((*active_)[ Gas ]);
|
||||
invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0;
|
||||
}
|
||||
|
||||
// Store the perforation phase flux for later usage.
|
||||
@ -253,7 +271,7 @@ namespace Opm {
|
||||
EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt;
|
||||
resWell_loc += getQs(w, p1);
|
||||
for (int p2 = 0; p2 < np; ++p2) {
|
||||
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+3);
|
||||
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+blocksize);
|
||||
}
|
||||
resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value();
|
||||
}
|
||||
@ -464,11 +482,11 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
typedef DenseAd::Evaluation<double, /*size=*/3> Eval;
|
||||
typedef DenseAd::Evaluation<double, /*size=*/blocksize> Eval;
|
||||
EvalWell extendEval(Eval in) const {
|
||||
EvalWell out = 0.0;
|
||||
out.setValue(in.value());
|
||||
for(int i = 0;i<3;++i) {
|
||||
for(int i = 0; i < blocksize;++i) {
|
||||
out.setDerivative(i, in.derivative(flowToEbosPvIdx(i)));
|
||||
}
|
||||
return out;
|
||||
@ -482,7 +500,7 @@ namespace Opm {
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
wellVariables_[w + nw*phaseIdx] = 0.0;
|
||||
wellVariables_[w + nw*phaseIdx].setValue(xw.wellSolutions()[w + nw* phaseIdx]);
|
||||
wellVariables_[w + nw*phaseIdx].setDerivative(np + phaseIdx, 1.0);
|
||||
wellVariables_[w + nw*phaseIdx].setDerivative(blocksize + phaseIdx, 1.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -819,10 +837,14 @@ namespace Opm {
|
||||
const PhaseUsage& pu = fluid_->phaseUsage();
|
||||
const int np = fluid_->numPhases();
|
||||
b_perf.resize(nperf*np);
|
||||
rsmax_perf.resize(nperf);
|
||||
rvmax_perf.resize(nperf);
|
||||
surf_dens_perf.resize(nperf*np);
|
||||
|
||||
//rs and rv are only used if both oil and gas is present
|
||||
if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) {
|
||||
rsmax_perf.resize(nperf);
|
||||
rvmax_perf.resize(nperf);
|
||||
}
|
||||
|
||||
// Compute the average pressure in each well block
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
|
||||
@ -915,34 +937,64 @@ namespace Opm {
|
||||
|
||||
// update the second and third well variable (The flux fractions)
|
||||
std::vector<double> F(np,0.0);
|
||||
const int sign2 = dwells[w][flowPhaseToEbosCompIdx(1)] > 0 ? 1: -1;
|
||||
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(1)]),dFLimit);
|
||||
well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited;
|
||||
const int sign3 = dwells[w][flowPhaseToEbosCompIdx(2)] > 0 ? 1: -1;
|
||||
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(2)]),dFLimit);
|
||||
well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited;
|
||||
F[Water] = well_state.wellSolutions()[nw + w];
|
||||
F[Gas] = well_state.wellSolutions()[2*nw + w];
|
||||
F[Oil] = 1.0 - F[Water] - F[Gas];
|
||||
|
||||
if (F[Water] < 0.0) {
|
||||
F[Gas] /= (1.0 - F[Water]);
|
||||
F[Oil] /= (1.0 - F[Water]);
|
||||
F[Water] = 0.0;
|
||||
if ((*active_)[ Water ]) {
|
||||
const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1;
|
||||
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit);
|
||||
well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited;
|
||||
}
|
||||
if (F[Gas] < 0.0) {
|
||||
F[Water] /= (1.0 - F[Gas]);
|
||||
F[Oil] /= (1.0 - F[Gas]);
|
||||
F[Gas] = 0.0;
|
||||
|
||||
if ((*active_)[ Gas ]) {
|
||||
const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1;
|
||||
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit);
|
||||
well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited;
|
||||
}
|
||||
|
||||
assert((*active_)[ Oil ]);
|
||||
F[Oil] = 1.0;
|
||||
if ((*active_)[ Water ]) {
|
||||
F[Water] = well_state.wellSolutions()[WFrac*nw + w];
|
||||
F[Oil] -= F[Water];
|
||||
}
|
||||
|
||||
if ((*active_)[ Gas ]) {
|
||||
F[Gas] = well_state.wellSolutions()[GFrac*nw + w];
|
||||
F[Oil] -= F[Gas];
|
||||
}
|
||||
|
||||
if ((*active_)[ Water ]) {
|
||||
if (F[Water] < 0.0) {
|
||||
if ((*active_)[ Gas ]) {
|
||||
F[Gas] /= (1.0 - F[Water]);
|
||||
}
|
||||
F[Oil] /= (1.0 - F[Water]);
|
||||
F[Water] = 0.0;
|
||||
}
|
||||
}
|
||||
if ((*active_)[ Gas ]) {
|
||||
if (F[Gas] < 0.0) {
|
||||
if ((*active_)[ Water ]) {
|
||||
F[Water] /= (1.0 - F[Gas]);
|
||||
}
|
||||
F[Oil] /= (1.0 - F[Gas]);
|
||||
F[Gas] = 0.0;
|
||||
}
|
||||
}
|
||||
if (F[Oil] < 0.0) {
|
||||
F[Water] /= (1.0 - F[Oil]);
|
||||
F[Gas] /= (1.0 - F[Oil]);
|
||||
if ((*active_)[ Water ]) {
|
||||
F[Water] /= (1.0 - F[Oil]);
|
||||
}
|
||||
if ((*active_)[ Gas ]) {
|
||||
F[Gas] /= (1.0 - F[Oil]);
|
||||
}
|
||||
F[Oil] = 0.0;
|
||||
}
|
||||
well_state.wellSolutions()[nw + w] = F[Water];
|
||||
well_state.wellSolutions()[2*nw + w] = F[Gas];
|
||||
//std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
|
||||
|
||||
if ((*active_)[ Water ]) {
|
||||
well_state.wellSolutions()[WFrac*nw + w] = F[Water];
|
||||
}
|
||||
if ((*active_)[ Gas ]) {
|
||||
well_state.wellSolutions()[GFrac*nw + w] = F[Gas];
|
||||
}
|
||||
|
||||
// The interpretation of the first well variable depends on the well control
|
||||
const WellControls* wc = wells().ctrls[w];
|
||||
@ -969,18 +1021,18 @@ namespace Opm {
|
||||
case THP: // The BHP and THP both uses the total rate as first well variable.
|
||||
case BHP:
|
||||
{
|
||||
well_state.wellSolutions()[w] = xvar_well_old[w] - dwells[w][flowPhaseToEbosCompIdx(0)];
|
||||
well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)];
|
||||
|
||||
switch (wells().type[w]) {
|
||||
case INJECTOR:
|
||||
for (int p = 0; p < np; ++p) {
|
||||
const double comp_frac = wells().comp_frac[np*w + p];
|
||||
well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w];
|
||||
well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[nw*XvarWell + w];
|
||||
}
|
||||
break;
|
||||
case PRODUCER:
|
||||
for (int p = 0; p < np; ++p) {
|
||||
well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p];
|
||||
well_state.wellRates()[w*np + p] = well_state.wellSolutions()[nw*XvarWell + w] * F[p];
|
||||
}
|
||||
break;
|
||||
}
|
||||
@ -1038,10 +1090,10 @@ namespace Opm {
|
||||
case SURFACE_RATE: // Both rate controls use bhp as first well variable
|
||||
case RESERVOIR_RATE:
|
||||
{
|
||||
const int sign1 = dwells[w][flowPhaseToEbosCompIdx(0)] > 0 ? 1: -1;
|
||||
const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(0)]),std::abs(xvar_well_old[w])*dBHPLimit);
|
||||
well_state.wellSolutions()[w] = std::max(xvar_well_old[w] - dx1_limited,1e5);
|
||||
well_state.bhp()[w] = well_state.wellSolutions()[w];
|
||||
const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1;
|
||||
const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit);
|
||||
well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5);
|
||||
well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w];
|
||||
|
||||
if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
|
||||
if (wells().type[w]==PRODUCER) {
|
||||
@ -1235,14 +1287,14 @@ namespace Opm {
|
||||
case BHP:
|
||||
{
|
||||
const WellType& well_type = wells().type[w];
|
||||
xw.wellSolutions()[w] = 0.0;
|
||||
xw.wellSolutions()[nw*XvarWell + w] = 0.0;
|
||||
if (well_type == INJECTOR) {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p];
|
||||
xw.wellSolutions()[nw*XvarWell + w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p];
|
||||
}
|
||||
} else {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p];
|
||||
xw.wellSolutions()[nw*XvarWell + w] += g[p] * xw.wellRates()[np*w + p];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1252,7 +1304,7 @@ namespace Opm {
|
||||
case RESERVOIR_RATE: // Intentional fall-through
|
||||
case SURFACE_RATE:
|
||||
{
|
||||
xw.wellSolutions()[w] = xw.bhp()[w];
|
||||
xw.wellSolutions()[nw*XvarWell + w] = xw.bhp()[w];
|
||||
}
|
||||
break;
|
||||
}
|
||||
@ -1262,11 +1314,20 @@ namespace Opm {
|
||||
tot_well_rate += g[p] * xw.wellRates()[np*w + p];
|
||||
}
|
||||
if(std::abs(tot_well_rate) > 0) {
|
||||
xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate;
|
||||
xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ;
|
||||
if ((*active_)[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate;
|
||||
}
|
||||
if ((*active_)[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ;
|
||||
}
|
||||
} else {
|
||||
xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water];
|
||||
xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas];
|
||||
if ((*active_)[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + w] = wells().comp_frac[np*w + Water];
|
||||
}
|
||||
|
||||
if ((*active_)[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + w] = wells().comp_frac[np*w + Gas];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1482,13 +1543,15 @@ namespace Opm {
|
||||
return bhp;
|
||||
|
||||
}
|
||||
return wellVariables_[wellIdx];
|
||||
const int nw = wells().number_of_wells;
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
|
||||
EvalWell getQs(const int wellIdx, const int phaseIdx) const {
|
||||
EvalWell qs = 0.0;
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
const int np = wells().number_of_phases;
|
||||
const int nw = wells().number_of_wells;
|
||||
const double target_rate = well_controls_get_current_target(wc);
|
||||
|
||||
if (wells().type[wellIdx] == INJECTOR) {
|
||||
@ -1497,7 +1560,7 @@ namespace Opm {
|
||||
return qs;
|
||||
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
|
||||
return wellVariables_[wellIdx];
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
qs.setValue(target_rate);
|
||||
return qs;
|
||||
@ -1505,7 +1568,7 @@ namespace Opm {
|
||||
|
||||
// Producers
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
|
||||
return wellVariables_[wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
}
|
||||
if (well_controls_get_current_type(wc) == SURFACE_RATE) {
|
||||
const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx];
|
||||
@ -1530,18 +1593,25 @@ namespace Opm {
|
||||
}
|
||||
|
||||
EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const {
|
||||
assert(wells().number_of_phases == 3);
|
||||
const int nw = wells().number_of_wells;
|
||||
if (phaseIdx == Water) {
|
||||
return wellVariables_[nw + wellIdx];
|
||||
return wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if (phaseIdx == Gas) {
|
||||
return wellVariables_[2*nw + wellIdx];
|
||||
return wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
// Oil
|
||||
return 1.0 - wellVariables_[nw + wellIdx] - wellVariables_[2 * nw + wellIdx];
|
||||
// Oil fraction
|
||||
EvalWell well_fraction = 1.0;
|
||||
if ((*active_)[Water]) {
|
||||
well_fraction -= wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if ((*active_)[Gas]) {
|
||||
well_fraction -= wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
return well_fraction;
|
||||
}
|
||||
|
||||
EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const {
|
||||
|
||||
@ -62,7 +62,7 @@ namespace Opm
|
||||
/// to give useful initial values to the bhp(), wellRates()
|
||||
/// and perfPhaseRates() fields, depending on controls
|
||||
template <class State, class PrevState>
|
||||
void init(const Wells* wells, const State& state, const PrevState& prevState)
|
||||
void init(const Wells* wells, const State& state, const PrevState& prevState, const PhaseUsage& pu)
|
||||
{
|
||||
// call init on base class
|
||||
BaseType :: init(wells, state, prevState);
|
||||
@ -115,19 +115,31 @@ namespace Opm
|
||||
}
|
||||
break;
|
||||
}
|
||||
#warning "HACK: that assert() was probably there for a reason!"
|
||||
//assert(np == 3);
|
||||
|
||||
double total_rates = 0.0;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
total_rates += g[p] * wellRates()[np*w + p];
|
||||
}
|
||||
|
||||
const int waterpos = pu.phase_pos[Water];
|
||||
const int gaspos = pu.phase_pos[Gas];
|
||||
|
||||
assert (np == 3 || np == 2 && !pu.phase_used[Gas] );
|
||||
// assumes the gas fractions are stored after water fractions
|
||||
if(std::abs(total_rates) > 0) {
|
||||
wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water;
|
||||
wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas
|
||||
if( pu.phase_used[Water] ) {
|
||||
wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates;
|
||||
}
|
||||
if( pu.phase_used[Gas] ) {
|
||||
wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + gaspos] / total_rates ;
|
||||
}
|
||||
} else {
|
||||
wellSolutions()[nw + w] = wells->comp_frac[np*w + Water];
|
||||
wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas];
|
||||
if( pu.phase_used[Water] ) {
|
||||
wellSolutions()[nw + w] = wells->comp_frac[np*w + waterpos];
|
||||
}
|
||||
if( pu.phase_used[Gas] ) {
|
||||
wellSolutions()[2*nw + w] = wells->comp_frac[np*w + gaspos];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -159,9 +171,9 @@ namespace Opm
|
||||
}
|
||||
|
||||
template <class State>
|
||||
void resize(const Wells* wells, const State& state) {
|
||||
void resize(const Wells* wells, const State& state, const PhaseUsage& pu ) {
|
||||
const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes
|
||||
init(wells, state, dummy_state) ;
|
||||
init(wells, state, dummy_state, pu) ;
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user