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 #1176 from totto82/impl_solvent

Browse Source 
Implement solvent model in flow_ebos
This commit is contained in:
Atgeirr Flø Rasmussen 2017-06-02 10:08:49 +02:00 committed by GitHub
commit 4c70c831d4
12 changed files with 376 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -48,7 +48,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/VFPInjProperties.cpp
opm/autodiff/WellMultiSegment.cpp
opm/autodiff/MultisegmentWells.cpp
opm/autodiff/BlackoilSolventState.cpp
opm/autodiff/MissingFeatures.cpp
opm/polymer/PolymerState.cpp
opm/polymer/PolymerBlackoilState.cpp
@ -169,7 +168,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/BlackoilSequentialModel.hpp
opm/autodiff/BlackoilSolventModel.hpp
opm/autodiff/BlackoilSolventModel_impl.hpp
opm/autodiff/BlackoilSolventState.hpp
opm/autodiff/BlackoilMultiSegmentModel.hpp
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
opm/autodiff/BlackoilTransportModel.hpp

99
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,8 +627,14 @@ 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);
maxVal = std::max(std::abs(dss),maxVal);
double step = dsMax()/maxVal;
step = std::min(step, 1.0);
@ -638,6 +647,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 +697,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 +885,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 +919,24 @@ 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();
R2_idx[cell_idx] = ebosResid[cell_idx][solventCompIdx];
}
}
Vector pv_vector;
@ -969,6 +1001,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 +1064,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 +1312,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 +1400,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 +1513,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 +1574,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 +1653,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 +1697,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 +1729,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 ];
}

3
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -22,7 +22,6 @@
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilSolventState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/StandardWellsSolvent.hpp>
@ -255,7 +254,7 @@ namespace Opm {
template <class Grid>
struct ModelTraits< BlackoilSolventModel<Grid> >
{
typedef BlackoilSolventState ReservoirState;
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoilSolvent WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilSolventSolutionState SolutionState;

2
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -135,7 +135,7 @@ namespace Opm {
// Initial solvent concentration.
if (has_solvent_) {
const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL );
const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL );
const int nc = solvent_saturation.size();
const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);

34
opm/autodiff/BlackoilSolventState.cpp
View File

@ -1,34 +0,0 @@
/*
Copyright 2015 IRIS AS, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/BlackoilSolventState.hpp>
namespace Opm
{
const std::string BlackoilSolventState::SSOL = "SSOL";
BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases)
: BlackoilState( number_of_cells , number_of_faces , number_of_phases)
{
registerCellData( SSOL , 1 );
}
} // namespace Opm

46
opm/autodiff/BlackoilSolventState.hpp
View File

@ -1,46 +0,0 @@
/*
Copyright 2015 IRIS AS, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
#define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
#include <string>
#include <opm/core/simulator/BlackoilState.hpp>
namespace Opm
{
/// Simulator state for blackoil simulator with solvent.
/// We use the Blackoil state parameters.
class BlackoilSolventState : public BlackoilState
{
public:
static const std::string SSOL;
BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases);
};
} // namespace Opm
#endif // OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED

3
opm/autodiff/Compat.cpp
View File

@ -29,7 +29,6 @@
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
#include <opm/autodiff/BlackoilSolventState.hpp>
#include <opm/output/data/Cells.hpp>
#include <opm/output/data/Solution.hpp>
#include <opm/output/data/Wells.hpp>
@ -111,7 +110,7 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir,
}
if (phases.has_solvent) {
sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilSolventState::SSOL ) , data::TargetType::RESTART_SOLUTION );
sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION );
}
return sol;

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());
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;

3
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
View File

@ -22,7 +22,6 @@
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/autodiff/BlackoilSolventState.hpp>
#include <opm/autodiff/BlackoilSolventModel.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
@ -85,7 +84,7 @@ namespace Opm
struct SimulatorTraits<SimulatorFullyImplicitBlackoilSolvent<GridT> >
{
typedef WellStateFullyImplicitBlackoilSolvent WellState;
typedef BlackoilSolventState ReservoirState;
typedef BlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilSolventModel<Grid> Model;

22
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,15 @@ 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 +288,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 +397,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;
};

204
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
@ -1129,14 +1157,14 @@ namespace Opm {
if (pu.phase_used[BlackoilPhases::Vapour]) {
const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp;
const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp;
const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
if (pu.phase_used[BlackoilPhases::Liquid]) {
const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp;
const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
if (oilrate > 0) {
const double gasrate = std::abs(xw.wellRates()[gaspos_well]);
const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
double rv = 0.0;
if (gasrate > 0) {
rv = oilrate / gasrate;
@ -1156,11 +1184,11 @@ namespace Opm {
if (pu.phase_used[BlackoilPhases::Liquid]) {
const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp;
const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp;
const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
if (pu.phase_used[BlackoilPhases::Vapour]) {
rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp;
const double gasrate = std::abs(xw.wellRates()[gaspos_well]);
const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
if (gasrate > 0) {
const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
double rs = 0.0;
@ -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());
}
// We 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,15 @@ namespace Opm {
if (active_[ Gas ]) {
well_state.wellSolutions()[GFrac*nw + w] = F[Gas];
}
if(has_solvent_) {
well_state.wellSolutions()[SFrac*nw + w] = F_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.
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 +1763,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 +2141,34 @@ 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_ ) {
double comp_frac = 0.0;
if (compIdx == solventCompIdx) { // solvent
comp_frac = wells().comp_frac[np*wellIdx + pu.phase_pos[ Gas ]] * wsolvent(wellIdx);
} else if (compIdx == pu.phase_pos[ Gas ]) {
comp_frac = wells().comp_frac[np*wellIdx + compIdx] * (1.0 - wsolvent(wellIdx));
} else {
comp_frac = wells().comp_frac[np*wellIdx + compIdx];
}
if (comp_frac == 0.0) {
return qs; //zero
}
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
return comp_frac * wellVariables_[nw*XvarWell + wellIdx];
}
qs.setValue(comp_frac * target_rate);
return qs;
}
const double comp_frac = wells().comp_frac[np*wellIdx + compIdx];
if (comp_frac == 0.0) {
return qs;
@ -2117,17 +2214,27 @@ namespace Opm {
assert(phase_under_control >= 0);
EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(wellIdx, phase_under_control);
if (has_solvent_ && phase_under_control == Gas) {
// for GRAT controlled wells solvent is included in the target
wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventCompIdx);
}
if (compIdx == phase_under_control) {
if (has_solvent_ && phase_under_control == Gas) {
qs.setValue(target_rate * wellVolumeFractionScaled(wellIdx, Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() );
return qs;
}
qs.setValue(target_rate);
return qs;
}
// TODO: not sure why the single phase under control will have near zero fraction
const double eps = 1e-6;
if (wellVolumeFractionScaled(wellIdx, phase_under_control) < eps) {
if (wellVolumeFractionScaledPhaseUnderControl < eps) {
return qs;
}
return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaled(wellIdx, phase_under_control));
return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaledPhaseUnderControl);
}
// when it is a combined two phase rate limit, such like LRAT
@ -2175,6 +2282,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 +2295,9 @@ namespace Opm {
if (active_[Gas]) {
well_fraction -= wellVariables_[GFrac * nw + wellIdx];
}
if (has_solvent_) {
well_fraction -= wellVariables_[SFrac * nw + wellIdx];
}
return well_fraction;
}
@ -2198,6 +2312,10 @@ namespace Opm {
{
const WellControls* wc = wells().ctrls[wellIdx];
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
if (has_solvent_ && compIdx == solventCompIdx) {
return wellVolumeFraction(wellIdx, compIdx);
}
const double* distr = well_controls_get_current_distr(wc);
assert(compIdx < 3);
if (distr[compIdx] > 0.) {
@ -2208,8 +2326,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 +2718,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 +2737,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;
}
@ -2672,6 +2795,7 @@ namespace Opm {
std::vector<double>& well_flux) const
{
const int np = wells().number_of_phases;
const int numComp = numComponents();
well_flux.resize(np, 0.0);
const bool allow_cf = allow_cross_flow(well_index, ebosSimulator);
@ -2679,8 +2803,8 @@ namespace Opm {
const int cell_index = wells().well_cells[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0));
// flux for each perforation
std::vector<EvalWell> cq_s(np, 0.0);
std::vector<EvalWell> mob(np, 0.0);
std::vector<EvalWell> cq_s(numComp, 0.0);
std::vector<EvalWell> mob(numComp, 0.0);
getMobility(ebosSimulator, perf, cell_index, mob);
computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp,
wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
@ -2879,4 +3003,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;
assert(false);
return 0.0;
}
} // namespace Opm

87
opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
View File

@ -67,6 +67,48 @@ namespace Opm
// call init on base class
BaseType :: init(wells, state, prevState);
const int nw = wells->number_of_wells;
if (nw == 0) {
return;
}
const int nperf = wells->well_connpos[nw];
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
if (pu.has_solvent) {
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if( ! prevState.wellMap().empty() )
{
typedef typename WellMapType :: const_iterator const_iterator;
const_iterator end = prevState.wellMap().end();
for (int w = 0; w < nw; ++w) {
std::string name( wells->name[ w ] );
const_iterator it = prevState.wellMap().find( name );
if( it != end )
{
const int newIndex = w;
// perfSolventRates
int oldPerf_idx = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ];
const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
if( num_perf_old_well == num_perf_this_well )
{
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
}
}
}
}
}
}
// TODO: the reason to keep this is to avoid getting defaulted value BHP
// some facilities needed from opm-parser or opm-core
// It is a little tricky, since sometimes before applying group control, the only
@ -88,9 +130,12 @@ namespace Opm
if (nw == 0) {
return;
}
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 +181,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 +196,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 +216,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