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

Implement solvent model in flow_ebos

Browse Source 
1) Extends the well model to account for solvent surface volumes
2) Add solvent to updateState
3) Add solvent to well and field output

The solvent parts is encapsled in if (has_solvent_) and should not effect
the standard runs.
This commit is contained in:
Tor Harald Sandve 2017-05-09 08:21:51 +02:00
parent 138eb2c91b
commit b987e4b324
6 changed files with 298 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

98
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -84,6 +84,7 @@ NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
SET_BOOL_PROP(EclFlowProblem, DisableWells, true);
SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false);
SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true);
SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
// SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy
// code for fluid and satfunc handling gets fully retired.
@ -119,6 +120,7 @@ namespace Opm {
typedef double Scalar;
static const int numEq = BlackoilIndices::numEq;
static const int solventCompIdx = 3; //TODO get this from ebos
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@ -162,6 +164,7 @@ namespace Opm {
, active_(detail::activePhases(fluid.phaseUsage()))
, has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil())
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, param_( param )
, well_model_ (well_model)
, terminal_output_ (terminal_output)
@ -624,6 +627,10 @@ namespace Opm {
}
dso -= dsg;
// solvent
const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0;
dso -= dss;
// Appleyard chop process.
maxVal = std::max(std::abs(dsg),maxVal);
double step = dsMax()/maxVal;
@ -638,6 +645,12 @@ namespace Opm {
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
sg -= step * dsg;
}
if (has_solvent_) {
double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
ss -= step * dss;
}
double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
so -= step * dso;
@ -682,12 +695,20 @@ namespace Opm {
if (sg <= 0.0 && has_disgas_) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
sg = 0;
so = 1.0 - sw - sg;
so = 1.0 - sw;
if (has_solvent_) {
double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
so -= ss;
}
rs *= (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;
sg = 1.0 - sw;
if (has_solvent_) {
double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx];
sg -= ss;
}
rv *= (1-epsilon);
}
break;
@ -862,6 +883,7 @@ namespace Opm {
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = numPhases();
const int numComp = numComponents();
Vector R_sum(numComp);
@ -895,16 +917,25 @@ namespace Opm {
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx )
{
Vector& R2_idx = R2[ compIdx ];
Vector& B_idx = B[ compIdx ];
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(compIdx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx);
Vector& R2_idx = R2[ phaseIdx ];
Vector& B_idx = B[ phaseIdx ];
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx);
B_idx [cell_idx] = 1.0 / fs.invB(ebosPhaseIdx).value();
R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx];
}
if (has_solvent_ ) {
Vector& R2_idx = R2[ solventCompIdx ];
Vector& B_idx = B[ solventCompIdx ];
B_idx [cell_idx] = 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
const int ebosCompIdx = flowPhaseToEbosCompIdx(solventCompIdx);
R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx];
}
}
Vector pv_vector;
@ -969,6 +1000,9 @@ namespace Opm {
const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
key[ phaseIdx ] = std::toupper( phaseName.front() );
}
if (has_solvent_) {
key[ solventCompIdx ] = "S";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " MB(" + key[ compIdx ] + ") ";
@ -1029,7 +1063,11 @@ namespace Opm {
if (numPhases() == 2) {
return 2;
}
return FluidSystem::numComponents;
int numComp = FluidSystem::numComponents;
if (has_solvent_)
numComp ++;
return numComp;
}
/// Wrapper required due to not following generic API
@ -1273,6 +1311,11 @@ namespace Opm {
VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" );
VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" );
if (has_solvent_) {
simData.registerCellData( "SSOL", 1 );
}
VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero;
std::vector<int> failed_cells_pb;
std::vector<int> failed_cells_pd;
const auto& gridView = ebosSimulator().gridView();
@ -1356,6 +1399,11 @@ namespace Opm {
krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
}
if (has_solvent_)
{
ssol[cellIdx] = intQuants.solventSaturation().value();
}
// hack to make the intial output of rs and rv Ecl compatible.
// For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
// where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
@ -1464,6 +1512,7 @@ namespace Opm {
const std::vector<int> cells_; // All grid cells
const bool has_disgas_;
const bool has_vapoil_;
const bool has_solvent_;
ModelParameters param_;
SimulatorReport failureReport_;
@ -1524,6 +1573,10 @@ namespace Opm {
// set water saturation
cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
if (has_solvent_) {
cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
}
// set switching variable and interpretation
if (active_[Gas] ) {
if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
@ -1599,19 +1652,18 @@ namespace Opm {
int flowToEbosPvIdx( const int flowPv ) const
{
assert(flowPv < 3);
const int flowToEbos[ 3 ] = {
const int flowToEbos[ 4 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx
BlackoilIndices::compositionSwitchIdx,
BlackoilIndices::solventSaturationIdx
};
return flowToEbos[ flowPv ];
}
int flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
assert(phaseIdx < 3);
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx };
const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
return phaseToComp[ phaseIdx ];
}
@ -1644,6 +1696,13 @@ namespace Opm {
cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] /= refDens;
cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] *= cellVolume;
}
if (has_solvent_) {
// no need to store refDens for all cells?
const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
const auto& refDens = intQuants->solventRefDensity();
cellRes[ solventCompIdx ] /= refDens;
cellRes[ solventCompIdx ] *= cellVolume;
}
}
for( auto row = ebosJac.begin(); row != endrow; ++row )
@ -1669,13 +1728,24 @@ namespace Opm {
(*col)[ebosCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
}
}
if (has_solvent_) {
// TODO store refDens pr pvtRegion?
const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(rowIdx, /*timeIdx=*/0);
const auto& refDens = intQuants->solventRefDensity();
for( int pvIdx=0; pvIdx < numEq; ++pvIdx )
{
(*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
(*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
}
}
}
}
}
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
assert(phaseIdx < 3);
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
return flowToEbos[ phaseIdx ];
}

4
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -252,7 +252,9 @@ public:
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
WellModel well_model(wells, &(wells_manager.wellCollection()), model_param_, terminal_output_);
const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum());
const WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum());
auto solver = createSolver(well_model);
std::vector<std::vector<double>> currentFluidInPlace;

2
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -62,7 +62,7 @@ namespace Opm
{
class SimulationDataContainer;
class BlackoilState;
class BlackoilSolventState;
void outputStateVtk(const UnstructuredGrid& grid,
const Opm::SimulationDataContainer& state,

21
opm/autodiff/StandardWellsDense.hpp
View File

@ -62,7 +62,8 @@ namespace Opm {
enum WellVariablePositions {
XvarWell = 0,
WFrac = 1,
GFrac = 2
GFrac = 2,
SFrac = 3
};
@ -85,6 +86,7 @@ enum WellVariablePositions {
typedef double Scalar;
static const int numEq = BlackoilIndices::numEq;
static const int numWellEq = numEq; //number of wellEq is the same as numEq in the model
static const int solventCompIdx = 3; //TODO get this from ebos
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@ -102,8 +104,10 @@ enum WellVariablePositions {
// --------- Public methods ---------
StandardWellsDense(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const bool terminal_output);
const bool terminal_output,
const int current_index);
void init(const PhaseUsage phase_usage_arg,
const std::vector<bool>& active_arg,
@ -118,10 +122,14 @@ enum WellVariablePositions {
/// The number of components in the model.
int numComponents() const
{
if (phase_usage_.num_phases == 2) {
if (numPhases() == 2) {
return 2;
}
return FluidSystem::numComponents;
int numComp = FluidSystem::numComponents;
if (has_solvent_)
numComp ++;
return numComp;
}
@ -279,12 +287,15 @@ enum WellVariablePositions {
protected:
bool wells_active_;
const Wells* wells_;
const std::vector< const Well* > wells_ecl_;
// Well collection is used to enforce the group control
WellCollection* well_collection_;
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
int current_timeIdx_;
PhaseUsage phase_usage_;
std::vector<bool> active_;
@ -385,6 +396,8 @@ enum WellVariablePositions {
const int well_index,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential) const;
double wsolvent(const int well_index) const;
};

161
opm/autodiff/StandardWellsDense_impl.hpp
View File

@ -7,13 +7,18 @@ namespace Opm {
StandardWellsDense<TypeTag>::
StandardWellsDense(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const bool terminal_output)
const bool terminal_output,
const int current_timeIdx)
: wells_active_(wells_arg!=nullptr)
, wells_(wells_arg)
, wells_ecl_(wells_ecl)
, well_collection_(well_collection)
, param_(param)
, terminal_output_(terminal_output)
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, current_timeIdx_(current_timeIdx)
, well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0)
, well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0)
, well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0)
@ -168,6 +173,7 @@ namespace Opm {
{
const int nw = wells().number_of_wells;
const int numComp = numComponents();
const int np = numPhases();
// clear all entries
duneB_ = 0.0;
@ -223,7 +229,11 @@ namespace Opm {
}
// Store the perforation phase flux for later usage.
well_state.perfPhaseRates()[perf*numComp + componentIdx] = cq_s[componentIdx].value();
if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value();
} else {
well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value();
}
}
// Store the perforation pressure for later usage.
@ -267,6 +277,9 @@ namespace Opm {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
}
if (has_solvent_) {
mob[solventCompIdx] = extendEval(intQuants.solventMobility());
}
} else {
const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
@ -281,6 +294,11 @@ namespace Opm {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
}
// this may not work if viscosity and relperms has been modified?
if (has_solvent_) {
OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
}
}
}
@ -443,8 +461,7 @@ namespace Opm {
StandardWellsDense<TypeTag>::
flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
assert(phaseIdx < 3);
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx };
const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
return phaseToComp[ phaseIdx ];
}
@ -457,11 +474,11 @@ namespace Opm {
StandardWellsDense<TypeTag>::
flowToEbosPvIdx( const int flowPv ) const
{
assert(flowPv < 3);
const int flowToEbos[ 3 ] = {
const int flowToEbos[ 4 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx
BlackoilIndices::compositionSwitchIdx,
BlackoilIndices::solventSaturationIdx
};
return flowToEbos[ flowPv ];
}
@ -754,6 +771,9 @@ namespace Opm {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
}
if (has_solvent_) {
b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
}
// Pressure drawdown (also used to determine direction of flow)
EvalWell well_pressure = bhp + cdp;
@ -803,6 +823,10 @@ namespace Opm {
volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos];
}
if (has_solvent_) {
volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx];
}
if (active_[Oil] && active_[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
@ -989,6 +1013,10 @@ namespace Opm {
B += 1 / fs.invB(ebosPhaseIdx).value();
}
if (has_solvent_) {
auto& B = B_avg[ solventCompIdx ];
B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
}
}
// compute global average
@ -1181,6 +1209,12 @@ namespace Opm {
for (int p = 0; p < pu.num_phases; ++p) {
surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex());
}
#warning HACK use cell values for solvent injector
if (has_solvent_) {
b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value();
surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity();
}
}
}
}
@ -1220,6 +1254,12 @@ namespace Opm {
well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited;
}
if (has_solvent_) {
const int sign4 = dwells[w][flowPhaseToEbosCompIdx(SFrac)] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(SFrac)]),dFLimit);
well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited;
}
assert(active_[ Oil ]);
F[Oil] = 1.0;
if (active_[ Water ]) {
@ -1232,11 +1272,20 @@ namespace Opm {
F[Oil] -= F[Gas];
}
double F_solvent = 0.0;
if (has_solvent_) {
F_solvent = well_state.wellSolutions()[SFrac*nw + w];
F[Oil] -= F_solvent;
}
if (active_[ Water ]) {
if (F[Water] < 0.0) {
if (active_[ Gas ]) {
F[Gas] /= (1.0 - F[Water]);
}
if (has_solvent_) {
F_solvent /= (1.0 - F[Water]);
}
F[Oil] /= (1.0 - F[Water]);
F[Water] = 0.0;
}
@ -1246,6 +1295,9 @@ namespace Opm {
if (active_[ Water ]) {
F[Water] /= (1.0 - F[Gas]);
}
if (has_solvent_) {
F_solvent /= (1.0 - F[Gas]);
}
F[Oil] /= (1.0 - F[Gas]);
F[Gas] = 0.0;
}
@ -1257,6 +1309,9 @@ namespace Opm {
if (active_[ Gas ]) {
F[Gas] /= (1.0 - F[Oil]);
}
if (has_solvent_) {
F_solvent /= (1.0 - F[Oil]);
}
F[Oil] = 0.0;
}
@ -1266,6 +1321,14 @@ namespace Opm {
if (active_[ Gas ]) {
well_state.wellSolutions()[GFrac*nw + w] = F[Gas];
}
if(has_solvent_) {
well_state.wellSolutions()[SFrac*nw + w] = F_solvent;
}
# warning 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 output, wellControls, well groups, THP etc.
if (has_solvent_){
F[Gas] += F_solvent;
}
// The interpretation of the first well variable depends on the well control
const WellControls* wc = wells().ctrls[w];
@ -1699,9 +1762,21 @@ namespace Opm {
const double grav)
{
// Compute densities
const int nperf = depth_perf.size();
const int numComponent = b_perf.size() / nperf;
const int np = wells().number_of_phases;
std::vector<double> perfRates(b_perf.size(),0.0);
for (int perf = 0; perf < nperf; ++perf) {
for (int phase = 0; phase < np; ++phase) {
perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase];
}
if(has_solvent_) {
perfRates[perf*numComponent + solventCompIdx] = xw.perfRateSolvent()[perf];
}
}
well_perforation_densities_ =
WellDensitySegmented::computeConnectionDensities(
wells(), phase_usage_, xw.perfPhaseRates(),
wells(), phase_usage_, perfRates,
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// Compute pressure deltas
@ -2065,13 +2140,27 @@ namespace Opm {
EvalWell qs = 0.0;
const WellControls* wc = wells().ctrls[wellIdx];
const int np = wells().number_of_phases;
assert(compIdx < np);
assert(compIdx < numComponents());
const int nw = wells().number_of_wells;
const auto pu = phase_usage_;
const double target_rate = well_controls_get_current_target(wc);
// TODO: the formulation for the injectors decides it only work with single phase
// surface rate injection control. Improvement will be required.
if (wells().type[wellIdx] == INJECTOR) {
if (has_solvent_ ) {
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
OPM_THROW(std::runtime_error,"BHP controlled solvent injector is unsupported. Check well "
<< wells().name [wellIdx] );
}
if (compIdx == pu.phase_pos[ Gas ]) { //gas
qs.setValue(target_rate * (1.0 - wsolvent(wellIdx)));
return qs;
} else if (compIdx == solventCompIdx) { // solvent
qs.setValue(wsolvent(wellIdx) * target_rate);
return qs;
}
}
const double comp_frac = wells().comp_frac[np*wellIdx + compIdx];
if (comp_frac == 0.0) {
return qs;
@ -2175,6 +2264,10 @@ namespace Opm {
return wellVariables_[GFrac * nw + wellIdx];
}
if (compIdx == solventCompIdx) {
return wellVariables_[SFrac * nw + wellIdx];
}
// Oil fraction
EvalWell well_fraction = 1.0;
if (active_[Water]) {
@ -2184,6 +2277,9 @@ namespace Opm {
if (active_[Gas]) {
well_fraction -= wellVariables_[GFrac * nw + wellIdx];
}
if (has_solvent_) {
well_fraction -= wellVariables_[SFrac * nw + wellIdx];
}
return well_fraction;
}
@ -2208,8 +2304,7 @@ namespace Opm {
return wellVolumeFraction(wellIdx, compIdx);
}
}
assert(compIdx < 3);
std::vector<double> g = {1,1,0.01};
std::vector<double> g = {1,1,0.01,0.01};
return (wellVolumeFraction(wellIdx, compIdx) / g[compIdx]);
}
@ -2601,7 +2696,10 @@ namespace Opm {
xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
}
if (active_[ Gas ]) {
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
}
if (has_solvent_) {
xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * wsolvent(well_index) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
}
} else {
const WellType& well_type = wells().type[well_index];
@ -2617,7 +2715,10 @@ namespace Opm {
if (active_[Gas]) {
if (distr[Gas] > 0.0) {
xw.wellSolutions()[GFrac * nw + well_index] = 1.0;
xw.wellSolutions()[GFrac * nw + well_index] = 1.0 - wsolvent(well_index);
if (has_solvent_) {
xw.wellSolutions()[SFrac * nw + well_index] = wsolvent(well_index);
}
} else {
xw.wellSolutions()[GFrac * nw + well_index] = 0.0;
}
@ -2879,4 +2980,38 @@ namespace Opm {
return potentials;
}
template<typename TypeTag>
double
StandardWellsDense<TypeTag>::
wsolvent(const int well_index) const {
if (!has_solvent_) {
return 0.0;
}
// loop over all wells until we find the well with the matching name
for (const auto& well : wells_ecl_) {
if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) {
continue;
}
WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_);
if (injection.injectorType == WellInjector::GAS) {
double solventFraction = well->getSolventFraction(current_timeIdx_);
// Look until we find the correct well
if (well->name() == wells().name[well_index]) {
return solventFraction;
}
}
}
// we didn't find it return 0;
// or should we throw?
return 0.0;
}
} // namespace Opm

47
opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
View File

@ -88,9 +88,14 @@ namespace Opm
if (nw == 0) {
return;
}
const int nperf = wells_->well_connpos[nw];
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
const int np = wells_->number_of_phases;
const int numComp = pu.has_solvent? np+1:np;
well_solutions_.clear();
well_solutions_.resize(nw * np, 0.0);
well_solutions_.resize(nw * numComp, 0.0);
std::vector<double> g = {1.0,1.0,0.01};
for (int w = 0; w < nw; ++w) {
WellControls* wc = wells_->ctrls[w];
@ -136,8 +141,12 @@ namespace Opm
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 ;
wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ;
}
if( pu.has_solvent) {
wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates;
}
} else {
@ -147,6 +156,10 @@ namespace Opm
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos];
}
if(pu.has_solvent) {
wellSolutions()[3*nw + w] = 0;
}
}
}
}
@ -163,14 +176,44 @@ namespace Opm
std::vector<double>& wellSolutions() { return well_solutions_; }
const std::vector<double>& wellSolutions() const { return well_solutions_; }
/// One rate pr well connection.
std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
/// One rate pr well
const double solventWellRate(const int w) const {
double solvent_well_rate = 0.0;
for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
solvent_well_rate += perfRateSolvent_[perf];
}
return solvent_well_rate;
}
data::Wells report(const PhaseUsage& pu) const override {
data::Wells res = BaseType::report(pu);
const int nw = WellState::numWells();
// If there are now wells numPhases throws a floating point
// exception.
if (nw == 0) {
return res;
}
if(pu.has_solvent) {
// add solvent component
for( int w = 0; w < nw; ++w ) {
using rt = data::Rates::opt;
res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
}
}
return res;
}
private:
std::vector<double> well_solutions_;
std::vector<double> perfRateSolvent_;
};
} // namespace Opm