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 #1263 from totto82/gasoil

Browse Source 
Make 2p gas oil compile and run
This commit is contained in:
Atgeirr Flø Rasmussen
2017-09-21 14:03:33 +02:00
committed by GitHub
parent c1024ce2d1 5839962e31
commit 83ca6d41e3
10 changed files with 180 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
compareECLFiles.cmake
View File

@@ -106,7 +106,6 @@ set(rel_tol 1e-5)
add_test_compareECLFiles(spe1 SPE1CASE2 flow ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1 SPE1CASE2 flow ${abs_tol} ${rel_tol} compareECLFiles "")
add_test_compareECLFiles(spe1 SPE1CASE2 flow_ebos ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1 SPE1CASE2 flow_ebos ${abs_tol} ${rel_tol} compareECLFiles "")
add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow ${abs_tol} ${rel_tol} compareECLFiles "" spe1) add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow ${abs_tol} ${rel_tol} compareECLFiles "" spe1)
add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_ebos ${abs_tol} ${rel_tol} compareECLFiles "" spe1)
add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_ebos_2p ${abs_tol} ${rel_tol} compareECLFiles "" spe1) add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_ebos_2p ${abs_tol} ${rel_tol} compareECLFiles "" spe1)
add_test_compareECLFiles(spe1 SPE1CASE2 flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1 SPE1CASE2 flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "")
add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "" spe1) add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "" spe1)

29
examples/flow.cpp
View File

@@ -52,12 +52,17 @@ namespace Properties {
// Twophase case // Twophase case
/////////////////////////////////// ///////////////////////////////////
NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem)); NEW_TYPE_TAG(EclFlowOilWaterProblem, INHERITS_FROM(EclFlowProblem));
//! The indices required by the model //! The indices required by the model
SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices, SET_TYPE_PROP(EclFlowOilWaterProblem, Indices,
Ewoms::BlackOilTwoPhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent)?1:0, GET_PROP_VALUE(TypeTag, EnablePolymer)?1:0, /*PVOffset=*/0>); Ewoms::BlackOilTwoPhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent)?1:0, GET_PROP_VALUE(TypeTag, EnablePolymer)?1:0, /*PVOffset=*/0, /*disabledCompIdx=*/2>);
NEW_TYPE_TAG(EclFlowGasOilProblem, INHERITS_FROM(EclFlowProblem));
//! The indices required by the model
SET_TYPE_PROP(EclFlowGasOilProblem, Indices,
Ewoms::BlackOilTwoPhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent)?1:0, GET_PROP_VALUE(TypeTag, EnablePolymer)?1:0, /*PVOffset=*/0, /*disabledCompIdx=*/1>);
/////////////////////////////////// ///////////////////////////////////
// Polymer case // Polymer case
/////////////////////////////////// ///////////////////////////////////
@@ -168,8 +173,22 @@ int main(int argc, char** argv)
// Twophase case // Twophase case
if( phases.size() == 2 ) { if( phases.size() == 2 ) {
Opm::FlowMainEbos<TTAG(EclFlowTwoPhaseProblem)> mainfunc; // oil-gas
return mainfunc.execute(argc, argv, deck, eclipseState ); if (phases.active( Opm::Phase::GAS ))
{
Opm::FlowMainEbos<TTAG(EclFlowGasOilProblem)> mainfunc;
return mainfunc.execute(argc, argv, deck, eclipseState );
}
// oil-water
else if ( phases.active( Opm::Phase::WATER ) )
{
Opm::FlowMainEbos<TTAG(EclFlowOilWaterProblem)> mainfunc;
return mainfunc.execute(argc, argv, deck, eclipseState );
}
else {
std::cerr << "No suitable configuration found, valid are Twophase (oilwater and oilgas), polymer, solvent, or blackoil" << std::endl;
return EXIT_FAILURE;
}
} }
// Polymer case // Polymer case
else if ( phases.active( Opm::Phase::POLYMER ) ) { else if ( phases.active( Opm::Phase::POLYMER ) ) {

2
examples/flow_ebos_2p.cpp
View File

@@ -38,7 +38,7 @@ namespace Properties {
NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem)); NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem));
//! The indices required by the model //! The indices required by the model
SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices, SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices,
Ewoms::BlackOilTwoPhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent)?1:0, GET_PROP_VALUE(TypeTag, EnablePolymer)?1:0, /*PVOffset=*/0>); Ewoms::BlackOilTwoPhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent)?1:0, GET_PROP_VALUE(TypeTag, EnablePolymer)?1:0, /*PVOffset=*/0, /*disabledCompIdx=*/2>);
}} }}
// ----------------- Main program ----------------- // ----------------- Main program -----------------

62
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -635,7 +635,7 @@ namespace Opm {
PrimaryVariables& priVars = solution[ cell_idx ]; PrimaryVariables& priVars = solution[ cell_idx ];
const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)]; const double& dp = dx[cell_idx][Indices::pressureSwitchIdx];
double& p = priVars[Indices::pressureSwitchIdx]; double& p = priVars[Indices::pressureSwitchIdx];
const double& dp_rel_max = dpMaxRel(); const double& dp_rel_max = dpMaxRel();
const int sign_dp = dp > 0 ? 1: -1; const int sign_dp = dp > 0 ? 1: -1;
@@ -643,9 +643,8 @@ namespace Opm {
p = std::max(p, 0.0); p = std::max(p, 0.0);
// Saturation updates. // Saturation updates.
const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0; const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
const int xvar_ind = active_[Water] ? 2 : 1; const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
const double dxvar = active_[Gas] ? dx[cell_idx][flowPhaseToEbosCompIdx(xvar_ind)] : 0.0;
double dso = 0.0; double dso = 0.0;
double dsg = 0.0; double dsg = 0.0;
@@ -1533,51 +1532,36 @@ namespace Opm {
public: public:
int ebosCompToFlowPhaseIdx( const int compIdx ) const
{
assert(compIdx < 3);
const int compToPhase[ 3 ] = { Oil, Water, Gas };
return compToPhase[ compIdx ];
}
int flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[] = {
Indices::pressureSwitchIdx,
Indices::waterSaturationIdx,
Indices::compositionSwitchIdx,
Indices::solventSaturationIdx
};
if (flowPv > 2 )
return flowPv;
return flowToEbos[ flowPv ];
}
int flowPhaseToEbosCompIdx( const int phaseIdx ) const int flowPhaseToEbosCompIdx( const int phaseIdx ) const
{ {
const int phaseToComp[] = { const auto& pu = phaseUsage_;
FluidSystem::waterCompIdx, if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
FluidSystem::oilCompIdx, return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
FluidSystem::gasCompIdx if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
}; return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (phaseIdx > 2 ) // for other phases return the index
return phaseIdx; return phaseIdx;
return phaseToComp[ phaseIdx ];
} }
private:
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{ {
const auto& pu = phaseUsage_;
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
return FluidSystem::waterPhaseIdx;
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
return FluidSystem::oilPhaseIdx;
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
return FluidSystem::gasPhaseIdx;
assert(phaseIdx < 3); assert(phaseIdx < 3);
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; // for other phases return the index
return flowToEbos[ phaseIdx ]; return phaseIdx;
} }
private:
void updateRateConverter() void updateRateConverter()
{ {

8
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -900,7 +900,9 @@ protected:
// set non-switching primary variables // set non-switching primary variables
PrimaryVariables& cellPv = solution[ cellIdx ]; PrimaryVariables& cellPv = solution[ cellIdx ];
// set water saturation // set water saturation
cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; if ( active[Water] ) {
cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
}
if (has_solvent) { if (has_solvent) {
cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx]; cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
@@ -944,7 +946,9 @@ protected:
/*storeViscosity=*/false, /*storeViscosity=*/false,
/*storeEnthalpy=*/false> SatOnlyFluidState; /*storeEnthalpy=*/false> SatOnlyFluidState;
SatOnlyFluidState fluidState; SatOnlyFluidState fluidState;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); if ( active[Water] ) {
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::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);

18
opm/autodiff/StandardWell.hpp
View File

@@ -47,20 +47,21 @@ namespace Opm
using typename Base::BlackoilIndices; using typename Base::BlackoilIndices;
using typename Base::PolymerModule; using typename Base::PolymerModule;
using Base::numEq;
// the positions of the primary variables for StandardWell // the positions of the primary variables for StandardWell
// there are three primary variables, the second and the third ones are F_w and F_g // there are three primary variables, the second and the third ones are F_w and F_g
// the first one can be total rate (G_t) or bhp, based on the control // the first one can be total rate (G_t) or bhp, based on the control
enum WellVariablePositions {
XvarWell = 0, static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0);
WFrac = 1, static const int XvarWell = 0;
GFrac = 2, static const int WFrac = gasoil? -1000: 1;
SFrac = 3 static const int GFrac = gasoil? 1: 2;
}; static const int SFrac = 3;
using typename Base::Scalar; using typename Base::Scalar;
using typename Base::ConvergenceReport; using typename Base::ConvergenceReport;
using Base::numEq;
using Base::has_solvent; using Base::has_solvent;
using Base::has_polymer; using Base::has_polymer;
@@ -158,7 +159,6 @@ namespace Opm
using Base::getAllowCrossFlow; using Base::getAllowCrossFlow;
using Base::phaseUsage; using Base::phaseUsage;
using Base::active; using Base::active;
using Base::flowToEbosPvIdx;
using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosPhaseIdx;
using Base::flowPhaseToEbosCompIdx; using Base::flowPhaseToEbosCompIdx;
using Base::numComponents; using Base::numComponents;
@@ -296,6 +296,8 @@ namespace Opm
void getMobility(const Simulator& ebosSimulator, void getMobility(const Simulator& ebosSimulator,
const int perf, const int perf,
std::vector<EvalWell>& mob) const; std::vector<EvalWell>& mob) const;
double scalingFactor(const int comp_idx) const;
}; };
} }

160
opm/autodiff/StandardWell_impl.hpp
View File

@@ -299,27 +299,13 @@ namespace Opm
StandardWell<TypeTag>:: StandardWell<TypeTag>::
wellVolumeFractionScaled(const int compIdx) const wellVolumeFractionScaled(const int compIdx) const
{ {
// TODO: we should be able to set the g for the well based on the control type
// instead of using explicit code for g all the times
const WellControls* wc = well_controls_;
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
if (has_solvent && compIdx == contiSolventEqIdx) { const double scal = scalingFactor(compIdx);
return wellVolumeFraction(compIdx); if (scal > 0)
} return wellVolumeFraction(compIdx) / scal;
const double* distr = well_controls_get_current_distr(wc);
assert(compIdx < 3);
if (distr[compIdx] > 0.) {
return wellVolumeFraction(compIdx) / distr[compIdx];
} else {
// TODO: not sure why return EvalWell(0.) causing problem here
// Probably due to the wrong Jacobians.
return wellVolumeFraction(compIdx);
}
}
std::vector<double> g = {1, 1, 0.01, 0.01}; // the scaling factor may be zero for RESV controlled wells.
return (wellVolumeFraction(compIdx) / g[compIdx]); return wellVolumeFraction(compIdx);
} }
@@ -331,11 +317,12 @@ namespace Opm
StandardWell<TypeTag>:: StandardWell<TypeTag>::
wellVolumeFraction(const int compIdx) const wellVolumeFraction(const int compIdx) const
{ {
if (compIdx == Water) { const auto pu = phaseUsage();
if (active()[Water] && compIdx == pu.phase_pos[Water]) {
return primary_variables_evaluation_[WFrac]; return primary_variables_evaluation_[WFrac];
} }
if (compIdx == Gas) { if (active()[Gas] && compIdx == pu.phase_pos[Gas]) {
return primary_variables_evaluation_[GFrac]; return primary_variables_evaluation_[GFrac];
} }
@@ -390,7 +377,7 @@ namespace Opm
EvalWell out = 0.0; EvalWell out = 0.0;
out.setValue(in.value()); out.setValue(in.value());
for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx))); out.setDerivative(eqIdx, in.derivative(eqIdx));
} }
return out; return out;
} }
@@ -585,17 +572,11 @@ namespace Opm
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
if (!only_wells) { if (!only_wells) {
// also need to consider the efficiency factor when manipulating the jacobians. // also need to consider the efficiency factor when manipulating the jacobians.
ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx);
duneB_[0][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
} }
} }
// add a trivial equation for the dummy phase for 2p cases (Only support water + oil)
if ( numComp < numWellEq ) {
assert(!active()[ Gas ]);
invDuneD_[0][0][Gas][Gas] = 1.0;
}
// Store the perforation phase flux for later usage. // Store the perforation phase flux for later usage.
if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
@@ -613,7 +594,7 @@ namespace Opm
} }
if (!only_wells) { if (!only_wells) {
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx); ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
} }
ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
} }
@@ -787,6 +768,7 @@ namespace Opm
const int np = number_of_phases_; const int np = number_of_phases_;
const double dBHPLimit = param.dbhp_max_rel_; const double dBHPLimit = param.dbhp_max_rel_;
const double dFLimit = param.dwell_fraction_max_; const double dFLimit = param.dwell_fraction_max_;
const auto pu = phaseUsage();
const std::vector<double> xvar_well_old = primary_variables_; const std::vector<double> xvar_well_old = primary_variables_;
@@ -811,68 +793,68 @@ namespace Opm
} }
assert(active()[ Oil ]); assert(active()[ Oil ]);
F[Oil] = 1.0; F[pu.phase_pos[Oil]] = 1.0;
if (active()[ Water ]) { if (active()[ Water ]) {
F[Water] = primary_variables_[WFrac]; F[pu.phase_pos[Water]] = primary_variables_[WFrac];
F[Oil] -= F[Water]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
} }
if (active()[ Gas ]) { if (active()[ Gas ]) {
F[Gas] = primary_variables_[GFrac]; F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
F[Oil] -= F[Gas]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
} }
double F_solvent = 0.0; double F_solvent = 0.0;
if (has_solvent) { if (has_solvent) {
F_solvent = primary_variables_[SFrac]; F_solvent = primary_variables_[SFrac];
F[Oil] -= F_solvent; F[pu.phase_pos[Oil]] -= F_solvent;
} }
if (active()[ Water ]) { if (active()[ Water ]) {
if (F[Water] < 0.0) { if (F[Water] < 0.0) {
if (active()[ Gas ]) { if (active()[ Gas ]) {
F[Gas] /= (1.0 - F[Water]); F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
} }
if (has_solvent) { if (has_solvent) {
F_solvent /= (1.0 - F[Water]); F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
} }
F[Oil] /= (1.0 - F[Water]); F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
F[Water] = 0.0; F[pu.phase_pos[Water]] = 0.0;
} }
} }
if (active()[ Gas ]) { if (active()[ Gas ]) {
if (F[Gas] < 0.0) { if (F[pu.phase_pos[Gas]] < 0.0) {
if (active()[ Water ]) { if (active()[ Water ]) {
F[Water] /= (1.0 - F[Gas]); F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
} }
if (has_solvent) { if (has_solvent) {
F_solvent /= (1.0 - F[Gas]); F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
} }
F[Oil] /= (1.0 - F[Gas]); F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
F[Gas] = 0.0; F[pu.phase_pos[Gas]] = 0.0;
} }
} }
if (F[Oil] < 0.0) { if (F[pu.phase_pos[Oil]] < 0.0) {
if (active()[ Water ]) { if (active()[ Water ]) {
F[Water] /= (1.0 - F[Oil]); F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
} }
if (active()[ Gas ]) { if (active()[ Gas ]) {
F[Gas] /= (1.0 - F[Oil]); F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
} }
if (has_solvent) { if (has_solvent) {
F_solvent /= (1.0 - F[Oil]); F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
} }
F[Oil] = 0.0; F[pu.phase_pos[Oil]] = 0.0;
} }
if (active()[ Water ]) { if (active()[ Water ]) {
primary_variables_[WFrac] = F[Water]; primary_variables_[WFrac] = F[pu.phase_pos[Water]];
} }
if (active()[ Gas ]) { if (active()[ Gas ]) {
primary_variables_[GFrac] = F[Gas]; primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
} }
if(has_solvent) { if(has_solvent) {
primary_variables_[SFrac] = F_solvent; primary_variables_[SFrac] = F_solvent;
@@ -881,7 +863,7 @@ namespace Opm
// F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. // 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. // More testing is needed to make sure this is correct for well groups and THP.
if (has_solvent){ if (has_solvent){
F[Gas] += F_solvent; F[pu.phase_pos[Gas]] += F_solvent;
} }
// The interpretation of the first well variable depends on the well control // The interpretation of the first well variable depends on the well control
@@ -892,19 +874,12 @@ namespace Opm
const int current = well_state.currentControls()[index_of_well_]; const int current = well_state.currentControls()[index_of_well_];
const double target_rate = well_controls_iget_target(wc, current); const double target_rate = well_controls_iget_target(wc, current);
std::vector<double> g = {1,1,0.01}; for (int p = 0; p < np; ++p) {
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { const double scal = scalingFactor(p);
const double* distr = well_controls_iget_distr(wc, current); if (scal > 0) {
for (int p = 0; p < np; ++p) { F[p] /= scal ;
if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value } else {
F[p] /= distr[p]; F[p] = 0.;
} else {
F[p] = 0.;
}
}
} else {
for (int p = 0; p < np; ++p) {
F[p] /= g[p];
} }
} }
@@ -1225,7 +1200,7 @@ namespace Opm
const int w = index_of_well_; const int w = index_of_well_;
//rs and rv are only used if both oil and gas is present //rs and rv are only used if both oil and gas is present
if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) {
rsmax_perf.resize(nperf); rsmax_perf.resize(nperf);
rvmax_perf.resize(nperf); rvmax_perf.resize(nperf);
} }
@@ -1874,13 +1849,7 @@ namespace Opm
const int well_index = index_of_well_; const int well_index = index_of_well_;
const WellControls* wc = well_controls_; const WellControls* wc = well_controls_;
const double* distr = well_controls_get_current_distr(wc); const double* distr = well_controls_get_current_distr(wc);
const auto pu = phaseUsage();
std::vector<double> g = {1.0, 1.0, 0.01};
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
for (int phase = 0; phase < np; ++phase) {
g[phase] = distr[phase];
}
}
switch (well_controls_get_current_type(wc)) { switch (well_controls_get_current_type(wc)) {
case THP: case THP:
@@ -1892,7 +1861,7 @@ namespace Opm
} }
} else { } else {
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
primary_variables_[XvarWell] += g[p] * well_state.wellRates()[np*well_index + p]; primary_variables_[XvarWell] += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
} }
} }
break; break;
@@ -1905,17 +1874,17 @@ namespace Opm
double tot_well_rate = 0.0; double tot_well_rate = 0.0;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
tot_well_rate += g[p] * well_state.wellRates()[np*well_index + p]; tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
} }
if(std::abs(tot_well_rate) > 0) { if(std::abs(tot_well_rate) > 0) {
if (active()[ Water ]) { if (active()[ Water ]) {
primary_variables_[WFrac] = g[Water] * well_state.wellRates()[np*well_index + Water] / tot_well_rate; primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
} }
if (active()[ Gas ]) { if (active()[ Gas ]) {
primary_variables_[GFrac] = g[Gas] * (well_state.wellRates()[np*well_index + Gas] - well_state.solventWellRate(well_index)) / tot_well_rate ; primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
} }
if (has_solvent) { if (has_solvent) {
primary_variables_[SFrac] = g[Gas] * well_state.solventWellRate(well_index) / tot_well_rate ; primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / tot_well_rate ;
} }
} else { // tot_well_rate == 0 } else { // tot_well_rate == 0
if (well_type_ == INJECTOR) { if (well_type_ == INJECTOR) {
@@ -1929,7 +1898,7 @@ namespace Opm
} }
if (active()[Gas]) { if (active()[Gas]) {
if (distr[Gas] > 0.0) { if (distr[pu.phase_pos[Gas]] > 0.0) {
primary_variables_[GFrac] = 1.0 - wsolvent(); primary_variables_[GFrac] = 1.0 - wsolvent();
if (has_solvent) { if (has_solvent) {
primary_variables_[SFrac] = wsolvent(); primary_variables_[SFrac] = wsolvent();
@@ -2053,4 +2022,33 @@ namespace Opm
return thp; return thp;
} }
template<typename TypeTag>
double
StandardWell<TypeTag>::scalingFactor(const int phaseIdx) const
{
const WellControls* wc = well_controls_;
const double* distr = well_controls_get_current_distr(wc);
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
if (has_solvent && phaseIdx == contiSolventEqIdx )
OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented");
return distr[phaseIdx];
}
const auto& pu = phaseUsage();
if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
return 1.0;
if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
return 1.0;
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
return 0.01;
if (has_solvent && phaseIdx == contiSolventEqIdx )
return 0.01;
// we should not come this far
assert(false);
return 1.0;
}
} }

18
opm/autodiff/StandardWellsDense_impl.hpp
View File

@@ -45,13 +45,13 @@ namespace Opm {
// has to be set always for the convergence check! // has to be set always for the convergence check!
global_nc_ = global_nc; global_nc_ = global_nc;
phase_usage_ = phase_usage_arg;
active_ = active_arg;
if ( ! localWellsActive() ) { if ( ! localWellsActive() ) {
return; return;
} }
phase_usage_ = phase_usage_arg;
active_ = active_arg;
calculateEfficiencyFactors(); calculateEfficiencyFactors();
#ifndef NDEBUG #ifndef NDEBUG
@@ -302,9 +302,17 @@ namespace Opm {
StandardWellsDense<TypeTag>:: StandardWellsDense<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{ {
const auto& pu = phase_usage_;
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
return FluidSystem::waterPhaseIdx;
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
return FluidSystem::oilPhaseIdx;
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
return FluidSystem::gasPhaseIdx;
assert(phaseIdx < 3); assert(phaseIdx < 3);
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; // for other phases return the index
return flowToEbos[ phaseIdx ]; return phaseIdx;
} }

2
opm/autodiff/WellInterface.hpp
View File

@@ -268,8 +268,6 @@ namespace Opm
int flowPhaseToEbosCompIdx( const int phaseIdx ) const; int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowToEbosPvIdx( const int flowPv ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
// TODO: it is dumplicated with StandardWellsDense // TODO: it is dumplicated with StandardWellsDense

54
opm/autodiff/WellInterface_impl.hpp
View File

@@ -224,47 +224,41 @@ namespace Opm
WellInterface<TypeTag>:: WellInterface<TypeTag>::
flowPhaseToEbosCompIdx( const int phaseIdx ) const flowPhaseToEbosCompIdx( const int phaseIdx ) const
{ {
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; const auto& pu = phaseUsage();
if (phaseIdx > 2 ) if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
return phaseIdx; return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
return phaseToComp[ phaseIdx ]; if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// for other phases return the index
return phaseIdx;
} }
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[ 3 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx
};
if (flowPv > 2 )
return flowPv;
return flowToEbos[ flowPv ];
}
template<typename TypeTag> template<typename TypeTag>
int int
WellInterface<TypeTag>:: WellInterface<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{ {
assert(phaseIdx < 3); const auto& pu = phaseUsage();
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; if (active()[Water] && pu.phase_pos[Water] == phaseIdx) {
return flowToEbos[ phaseIdx ]; return FluidSystem::waterPhaseIdx;
} }
if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) {
return FluidSystem::oilPhaseIdx;
}
if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) {
return FluidSystem::gasPhaseIdx;
}
assert(phaseIdx < 3);
// for other phases return the index
return phaseIdx;
}