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 #1220 from totto82/impl_polymer_simpleWellModel

Browse Source 
Add flow_ebos_polymer and flow_ebos_solvent
This commit is contained in:
Atgeirr Flø Rasmussen 2017-06-27 10:51:42 +02:00 committed by GitHub
commit bf0937fae1
11 changed files with 533 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -112,6 +112,8 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/flow_legacy.cpp
examples/flow_sequential.cpp
examples/flow_ebos.cpp
examples/flow_ebos_solvent.cpp
examples/flow_ebos_polymer.cpp
examples/flow_multisegment.cpp
examples/flow_solvent.cpp
examples/sim_2p_incomp.cpp
@ -132,6 +134,8 @@ list (APPEND PROGRAM_SOURCE_FILES
examples/sim_2p_incomp_ad.cpp
examples/sim_2p_comp_reorder.cpp
examples/flow_ebos.cpp
examples/flow_ebos_solvent.cpp
examples/flow_ebos_polymer.cpp
examples/flow_legacy.cpp
examples/flow_sequential.cpp
examples/flow_solvent.cpp

2
examples/flow_ebos.cpp
View File

@ -34,6 +34,6 @@
// ----------------- Main program -----------------
int main(int argc, char** argv)
{
Opm::FlowMainEbos mainfunc;
Opm::FlowMainEbos<TTAG(EclFlowProblem)> mainfunc;
return mainfunc.execute(argc, argv);
}

42
examples/flow_ebos_polymer.cpp Normal file
View File

@ -0,0 +1,42 @@
/*
Copyright 2017 IRIS AS
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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/material/densead/Evaluation.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <dune/grid/CpGrid.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/autodiff/FlowMainEbos.hpp>
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowPolymerProblem, INHERITS_FROM(EclFlowProblem));
SET_BOOL_PROP(EclFlowPolymerProblem, EnablePolymer, true);
}}
// ----------------- Main program -----------------
int main(int argc, char** argv)
{
Opm::FlowMainEbos<TTAG(EclFlowPolymerProblem)> mainfunc;
return mainfunc.execute(argc, argv);
}

42
examples/flow_ebos_solvent.cpp Normal file
View File

@ -0,0 +1,42 @@
/*
Copyright 2017 IRIS AS
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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/material/densead/Evaluation.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <dune/grid/CpGrid.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/autodiff/FlowMainEbos.hpp>
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowSolventProblem, INHERITS_FROM(EclFlowProblem));
SET_BOOL_PROP(EclFlowSolventProblem, EnableSolvent, true);
}}
// ----------------- Main program -----------------
int main(int argc, char** argv)
{
Opm::FlowMainEbos<TTAG(EclFlowSolventProblem)> mainfunc;
return mainfunc.execute(argc, argv);
}

95
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -83,7 +83,6 @@ 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.
@ -97,6 +96,7 @@ namespace Opm {
/// where gas can be dissolved in oil and vice versa. It
/// uses an industry-standard TPFA discretization with per-phase
/// upwind weighting of mobilities.
template <class TypeTag>
class BlackoilModelEbos
{
public:
@ -105,7 +105,6 @@ namespace Opm {
typedef WellStateFullyImplicitBlackoilDense WellState;
typedef BlackoilModelParameters ModelParameters;
typedef typename TTAG(EclFlowProblem) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
@ -118,7 +117,11 @@ namespace Opm {
typedef double Scalar;
static const int numEq = BlackoilIndices::numEq;
static const int solventCompIdx = 3; //TODO get this from ebos
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@ -149,7 +152,8 @@ namespace Opm {
const ModelParameters& param,
const StandardWellsDense<TypeTag>& well_model,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output)
const bool terminal_output
)
: ebosSimulator_(ebosSimulator)
, grid_(ebosSimulator_.gridManager().grid())
, istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
@ -161,8 +165,9 @@ namespace Opm {
, has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil())
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
, param_( param )
, well_model_ (well_model)
, well_model_ (well_model)
, terminal_output_ (terminal_output)
, rate_converter_(phaseUsage_, ebosSimulator_.problem().pvtRegionArray().empty()?nullptr:ebosSimulator_.problem().pvtRegionArray().data(), AutoDiffGrid::numCells(grid_), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
, current_relaxation_(1.0)
@ -615,6 +620,9 @@ namespace Opm {
const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0;
dso -= dss;
// polymer
const double dc = has_polymer_ ? dx[cell_idx][BlackoilIndices::polymerConcentrationIdx] : 0.0;
// Appleyard chop process.
maxVal = std::max(std::abs(dsg),maxVal);
maxVal = std::max(std::abs(dss),maxVal);
@ -637,6 +645,12 @@ namespace Opm {
ss -= step * dss;
}
if (has_polymer_) {
double& c = reservoir_state.getCellData( reservoir_state.POLYMER )[cell_idx];
c -= step * dc;
c = std::max(c, 0.0);
}
double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
so -= step * dso;
@ -821,7 +835,6 @@ namespace Opm {
const double tol_cnv = param_.tolerance_cnv_;
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = numPhases();
const int numComp = numComponents();
@ -867,10 +880,16 @@ namespace Opm {
}
if ( has_solvent_ ) {
B_avg[ solventCompIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
const auto R2 = ebosResid[cell_idx][solventCompIdx];
R_sum[ solventCompIdx ] += R2;
maxCoeff[ solventCompIdx ] = std::max( maxCoeff[ solventCompIdx ], std::abs( R2 ) / pvValue );
B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value();
const auto R2 = ebosResid[cell_idx][contiSolventEqIdx];
R_sum[ contiSolventEqIdx ] += R2;
maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue );
}
if (has_polymer_ ) {
B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value();
const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx];
R_sum[ contiPolymerEqIdx ] += R2;
maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
}
}
@ -938,7 +957,11 @@ namespace Opm {
key[ phaseIdx ] = std::toupper( phaseName.front() );
}
if (has_solvent_) {
key[ solventCompIdx ] = "S";
key[ solventSaturationIdx ] = "S";
}
if (has_polymer_) {
key[ polymerConcentrationIdx ] = "P";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
@ -1004,6 +1027,9 @@ namespace Opm {
if (has_solvent_)
numComp ++;
if (has_polymer_)
numComp ++;
return numComp;
}
@ -1253,6 +1279,11 @@ namespace Opm {
}
VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero;
if (has_polymer_) {
simData.registerCellData( "POLYMER", 1 );
}
VectorType& cpolymer = has_polymer_ ? simData.getCellData( "POLYMER" ) : zero;
std::vector<int> failed_cells_pb;
std::vector<int> failed_cells_pd;
const auto& gridView = ebosSimulator().gridView();
@ -1341,6 +1372,11 @@ namespace Opm {
ssol[cellIdx] = intQuants.solventSaturation().value();
}
if (has_polymer_)
{
cpolymer[cellIdx] = intQuants.polymerConcentration().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
@ -1450,6 +1486,7 @@ namespace Opm {
const bool has_disgas_;
const bool has_vapoil_;
const bool has_solvent_;
const bool has_polymer_;
ModelParameters param_;
SimulatorReport failureReport_;
@ -1513,6 +1550,11 @@ namespace Opm {
cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
}
if (has_polymer_) {
cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
}
// set switching variable and interpretation
if (active_[Gas] ) {
if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
@ -1588,18 +1630,23 @@ namespace Opm {
int flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[ 4 ] = {
const int flowToEbos[ 3 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx,
BlackoilIndices::solventSaturationIdx
BlackoilIndices::compositionSwitchIdx
};
if (flowPv > 2 )
return flowPv;
return flowToEbos[ flowPv ];
}
int flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
if (phaseIdx > 2 )
return phaseIdx;
return phaseToComp[ phaseIdx ];
}
@ -1636,8 +1683,11 @@ namespace Opm {
// 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;
cellRes[ contiSolventEqIdx ] /= refDens;
cellRes[ contiSolventEqIdx ] *= cellVolume;
}
if (has_polymer_) {
cellRes[ contiPolymerEqIdx ] *= cellVolume;
}
}
@ -1670,10 +1720,17 @@ namespace Opm {
const auto& refDens = intQuants->solventRefDensity();
for( int pvIdx=0; pvIdx < numEq; ++pvIdx )
{
(*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
(*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
(*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] /= refDens;
(*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
}
}
if (has_polymer_) {
for( int pvIdx=0; pvIdx < numEq; ++pvIdx )
{
(*col)[contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume;
}
}
}
}
}

11
opm/autodiff/Compat.cpp
View File

@ -23,7 +23,9 @@
#include <algorithm>
#include <cassert>
#include <iostream>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
@ -113,6 +115,15 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir,
sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION );
}
if (phases.has_polymer) {
if (reservoir.hasCellData( PolymerBlackoilState::CONCENTRATION )) { // compatibility with legacy polymer
sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( PolymerBlackoilState::CONCENTRATION ) , data::TargetType::RESTART_SOLUTION );
} else {
sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::POLYMER ) , data::TargetType::RESTART_SOLUTION );
}
}
return sol;
}

27
opm/autodiff/FlowMainEbos.hpp
View File

@ -61,6 +61,7 @@
namespace Opm
{
// The FlowMain class is the ebos based black-oil simulator.
template <class TypeTag>
class FlowMainEbos
{
enum FileOutputValue{
@ -73,7 +74,6 @@ namespace Opm
};
public:
typedef TTAG(EclFlowProblem) TypeTag;
typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator;
typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
@ -82,7 +82,7 @@ namespace Opm
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator;
typedef Opm::SimulatorFullyImplicitBlackoilEbos<TypeTag> Simulator;
typedef typename Simulator::ReservoirState ReservoirState;
typedef typename Simulator::OutputWriter OutputWriter;
@ -534,6 +534,27 @@ namespace Opm
}
initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL"));
// Get initial polymer concentration from ebos
if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
auto& cpolymer = state_->getCellData( state_->POLYMER );
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
cpolymer[c] = ebosProblem().polymerConcentration(c);
}
}
// Get initial solvent saturation from ebos
if (GET_PROP_VALUE(TypeTag, EnableSolvent)) {
auto& solvent = state_->getCellData( state_->SSOL );
auto& sat = state_->saturation();
const int np = props.numPhases();
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
solvent[c] = ebosProblem().solventSaturation(c);
sat[c * np + pu.phase_pos[Water]];
}
}
}
// Extract messages from parser.
@ -667,7 +688,7 @@ namespace Opm
// fis_solver_
void setupLinearSolver()
{
typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType;
typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
extractParallelGridInformationToISTL(grid(), parallel_information_);
fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );

4
opm/autodiff/FlowMainPolymer.hpp
View File

@ -89,10 +89,6 @@ namespace Opm
}
}
// Setup linear solver.
// Writes to:
// fis_solver_

22
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -42,14 +42,12 @@
namespace Opm {
class SimulatorFullyImplicitBlackoilEbos;
//class StandardWellsDense<FluidSystem>;
/// a simulator for the blackoil model
template<class TypeTag>
class SimulatorFullyImplicitBlackoilEbos
{
public:
typedef typename TTAG(EclFlowProblem) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@ -58,11 +56,12 @@ public:
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
typedef WellStateFullyImplicitBlackoilDense WellState;
typedef BlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef BlackoilModelEbos Model;
typedef BlackoilModelEbos<TypeTag> Model;
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolver<Model> Solver;
typedef StandardWellsDense<TypeTag> WellModel;
@ -255,7 +254,8 @@ public:
solver_timer.start();
const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum());
WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum());
WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_,
timer.currentStepNum());
auto solver = createSolver(well_model);
@ -413,11 +413,12 @@ public:
{ return ebosSimulator_.gridManager().grid(); }
protected:
void handleAdditionalWellInflow(SimulatorTimer& /* timer */,
void handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& /* wells_manager */,
WellState& /* well_state */,
const Wells* /* wells */)
{ }
{
}
std::unique_ptr<Solver> createSolver(WellModel& well_model)
{
@ -439,7 +440,8 @@ protected:
legacyDepth_,
legacyPoreVolume_,
rateConverter_.get(),
globalNumCells);
globalNumCells,
grid());
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
model_param_,
well_model,
@ -852,6 +854,10 @@ protected:
}
}
// Data.
Simulator& ebosSimulator_;

32
opm/autodiff/StandardWellsDense.hpp
View File

@ -57,6 +57,8 @@
#include <opm/simulators/WellSwitchingLogger.hpp>
#include <math.h>
namespace Opm {
enum WellVariablePositions {
@ -85,8 +87,11 @@ 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
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? 3:numEq; // //numEq; //number of wellEq is only for 3 for polymer
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@ -94,8 +99,7 @@ enum WellVariablePositions {
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
// For the conversion between the surface volume rate and resrevoir voidage rate
using RateConverterType = RateConverter::
@ -116,7 +120,8 @@ enum WellVariablePositions {
const std::vector<double>& depth_arg,
const std::vector<double>& pv_arg,
const RateConverterType* rate_converter,
long int global_nc);
long int global_nc,
const auto& grid);
/// The number of components in the model.
@ -146,6 +151,7 @@ enum WellVariablePositions {
void
getMobility(const Simulator& ebosSimulator,
const int w,
const int perf,
const int cell_idx,
std::vector<EvalWell>& mob) const;
@ -294,6 +300,7 @@ enum WellVariablePositions {
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;
int current_timeIdx_;
PhaseUsage phase_usage_;
@ -314,6 +321,13 @@ enum WellVariablePositions {
std::vector<double> well_perforation_densities_;
std::vector<double> well_perforation_pressure_diffs_;
std::vector<double> wpolymer_;
std::vector<double> wsolvent_;
std::vector<double> wells_rep_radius_;
std::vector<double> wells_perf_length_;
std::vector<double> wells_bore_diameter_;
std::vector<EvalWell> wellVariables_;
std::vector<double> F0_;
@ -397,6 +411,14 @@ enum WellVariablePositions {
const std::vector<double>& initial_potential) const;
double wsolvent(const int well_index) const;
double wpolymer(const int well_index) const;
void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
void computeRepRadiusPerfLength(const auto& grid);
};

339
opm/autodiff/StandardWellsDense_impl.hpp
View File

@ -18,6 +18,7 @@ namespace Opm {
, param_(param)
, terminal_output_(terminal_output)
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
, 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)
@ -47,7 +48,8 @@ namespace Opm {
const std::vector<double>& depth_arg,
const std::vector<double>& pv_arg,
const RateConverterType* rate_converter,
long int global_nc)
long int global_nc,
const auto& grid)
{
// has to be set always for the convergence check!
global_nc_ = global_nc;
@ -118,6 +120,14 @@ namespace Opm {
// resize temporary class variables
Cx_.resize( duneC_.N() );
invDrw_.resize( invDuneD_.N() );
if (has_polymer_)
{
if (PolymerModule::hasPlyshlog()) {
computeRepRadiusPerfLength(grid);
}
}
}
@ -194,7 +204,7 @@ namespace Opm {
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
std::vector<EvalWell> cq_s(numComp,0.0);
std::vector<EvalWell> mob(numComp, 0.0);
getMobility(ebosSimulator, perf, cell_idx, mob);
getMobility(ebosSimulator, w, perf, cell_idx, mob);
computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
@ -215,13 +225,19 @@ namespace Opm {
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
if (!only_wells) {
// 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);
duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
}
invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq);
}
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
if (!only_wells) {
// 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);
duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
}
}
// add trivial equation for 2p cases (Only support water + oil)
if (numComp == 2) {
assert(!active_[ Gas ]);
@ -229,13 +245,28 @@ namespace Opm {
}
// Store the perforation phase flux for later usage.
if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
if (has_solvent_ && componentIdx == solventSaturationIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value();
} else {
well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value();
}
}
if (has_polymer_) {
EvalWell cq_s_poly = cq_s[Water];
if (wells().type[w] == INJECTOR) {
cq_s_poly *= wpolymer(w);
} else {
cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
}
if (!only_wells) {
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx);
}
ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
}
}
// Store the perforation pressure for later usage.
well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf];
}
@ -249,8 +280,15 @@ namespace Opm {
}
resWell_[w][componentIdx] += resWell_loc.value();
}
// add trivial equation for polymer
if (has_polymer_) {
invDuneD_[w][w][contiPolymerEqIdx][polymerConcentrationIdx] = 1.0; //
}
}
// do the local inversion of D.
localInvert( invDuneD_ );
}
@ -259,7 +297,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag >::
getMobility(const Simulator& ebosSimulator, const int perf, const int cell_idx, std::vector<EvalWell>& mob) const
getMobility(const Simulator& ebosSimulator, const int w, const int perf, const int cell_idx, std::vector<EvalWell>& mob) const
{
const int np = wells().number_of_phases;
@ -278,7 +316,7 @@ namespace Opm {
mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
}
if (has_solvent_) {
mob[solventCompIdx] = extendEval(intQuants.solventMobility());
mob[solventSaturationIdx] = extendEval(intQuants.solventMobility());
}
} else {
@ -300,6 +338,51 @@ namespace Opm {
OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
}
}
// modify the water mobility if polymer is present
if (has_polymer_) {
// assume fully mixture for wells.
EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
if (wells().type[w] == INJECTOR) {
const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) );
}
if (PolymerModule::hasPlyshlog()) {
// compute the well water velocity with out shear effects.
const int numComp = numComponents();
bool allow_cf = allow_cross_flow(w, ebosSimulator);
const EvalWell& bhp = getBhp(w);
std::vector<EvalWell> cq_s(numComp,0.0);
computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf];
const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
const auto& scaledDrainageInfo =
materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
const Scalar& Swcr = scaledDrainageInfo.Swcr;
const EvalWell poro = extendEval(intQuants.porosity());
const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
// guard against zero porosity and no water
const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
if (PolymerModule::hasShrate()) {
// TODO Use the same conversion as for the reservoar equations.
// Need the "permeability" of the well?
// For now use the same formula as in legacy.
waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf];
}
EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
intQuants.pvtRegionIndex(),
waterVelocity);
// modify the mobility with the shear factor and recompute the well fluxes.
mob[ Water ] /= shearFactor;
}
}
}
@ -455,13 +538,31 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<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>
int
StandardWellsDense<TypeTag>::
flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx };
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
if (phaseIdx > 2 )
return phaseIdx;
return phaseToComp[ phaseIdx ];
}
@ -469,24 +570,6 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[ 4 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx,
BlackoilIndices::solventSaturationIdx
};
return flowToEbos[ flowPv ];
}
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
@ -497,10 +580,6 @@ namespace Opm {
return flowToEbos[ phaseIdx ];
}
template<typename TypeTag>
std::vector<double>
StandardWellsDense<TypeTag>::
@ -686,9 +765,9 @@ namespace Opm {
setWellVariables(const WellState& xw)
{
const int nw = wells().number_of_wells;
// for two-phase numComp < numEq
// for two-phase numComp < numWellEq
const int numComp = numComponents();
for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) {
for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) {
for (int w = 0; w < nw; ++w) {
const unsigned int idx = nw * eqIdx + w;
assert( idx < wellVariables_.size() );
@ -697,7 +776,7 @@ namespace Opm {
eval = 0.0;
eval.setValue( xw.wellSolutions()[ idx ] );
eval.setDerivative(numWellEq + eqIdx, 1.0);
eval.setDerivative(numEq + eqIdx, 1.0);
}
}
}
@ -765,7 +844,7 @@ namespace Opm {
b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
}
if (has_solvent_) {
b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
b_perfcells_dense[solventSaturationIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
}
// Pressure drawdown (also used to determine direction of flow)
@ -817,7 +896,7 @@ namespace Opm {
}
if (has_solvent_) {
volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx];
volumeRatio += cmix_s[solventSaturationIdx] / b_perfcells_dense[solventSaturationIdx];
}
if (active_[Oil] && active_[Gas]) {
@ -952,8 +1031,9 @@ namespace Opm {
const int nw = wells().number_of_wells;
const int numComp = numComponents();
std::vector<double> res(numComp*nw);
std::vector<double> res(numEq*nw, 0.0);
for( int compIdx = 0; compIdx < numComp; ++compIdx) {
for (int wellIdx = 0; wellIdx < nw; ++wellIdx) {
int idx = wellIdx + nw*compIdx;
res[idx] = resWell_[ wellIdx ][ compIdx ];
@ -1006,7 +1086,7 @@ namespace Opm {
B += 1 / fs.invB(ebosPhaseIdx).value();
}
if (has_solvent_) {
auto& B = B_avg[ solventCompIdx ];
auto& B = B_avg[ solventSaturationIdx ];
B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
}
}
@ -1204,8 +1284,8 @@ namespace Opm {
// 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();
b_perf[numComp*perf + solventSaturationIdx] = intQuants.solventInverseFormationVolumeFactor().value();
surf_dens_perf[numComp*perf + solventSaturationIdx] = intQuants.solventRefDensity();
}
}
}
@ -1764,7 +1844,7 @@ namespace Opm {
perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase];
}
if(has_solvent_) {
perfRates[perf*numComponent + solventCompIdx] = xw.perfRateSolvent()[perf];
perfRates[perf*numComponent + solventSaturationIdx] = xw.perfRateSolvent()[perf];
}
}
well_perforation_densities_ =
@ -2143,7 +2223,7 @@ namespace Opm {
if (wells().type[wellIdx] == INJECTOR) {
if (has_solvent_ ) {
double comp_frac = 0.0;
if (compIdx == solventCompIdx) { // solvent
if (has_solvent_ && compIdx == solventSaturationIdx) { // 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));
@ -2209,7 +2289,7 @@ namespace Opm {
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);
wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventSaturationIdx);
}
if (compIdx == phase_under_control) {
@ -2274,7 +2354,7 @@ namespace Opm {
return wellVariables_[GFrac * nw + wellIdx];
}
if (compIdx == solventCompIdx) {
if (has_solvent_ && compIdx == solventSaturationIdx) {
return wellVariables_[SFrac * nw + wellIdx];
}
@ -2305,7 +2385,7 @@ namespace Opm {
const WellControls* wc = wells().ctrls[wellIdx];
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
if (has_solvent_ && compIdx == solventCompIdx) {
if (has_solvent_ && compIdx == solventSaturationIdx) {
return wellVolumeFraction(wellIdx, compIdx);
}
const double* distr = well_controls_get_current_distr(wc);
@ -2710,10 +2790,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] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (xw.wellRates()[np*well_index + Gas] - xw.solventWellRate(well_index)) / 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 ;
xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * xw.solventWellRate(well_index) / tot_well_rate ;
}
} else {
const WellType& well_type = wells().type[well_index];
@ -2797,7 +2877,7 @@ namespace Opm {
// flux for each perforation
std::vector<EvalWell> cq_s(numComp, 0.0);
std::vector<EvalWell> mob(numComp, 0.0);
getMobility(ebosSimulator, perf, cell_index, mob);
getMobility(ebosSimulator, well_index, perf, cell_index, mob);
computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp,
wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
@ -3027,6 +3107,171 @@ namespace Opm {
}
template<typename TypeTag>
double
StandardWellsDense<TypeTag>::
wpolymer(const int well_index) const {
if (!has_polymer_) {
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_);
WellPolymerProperties polymer = well->getPolymerProperties(current_timeIdx_);
if (injection.injectorType == WellInjector::WATER) {
double polymerFraction = polymer.m_polymerConcentration;
// Look until we find the correct well
if (well->name() == wells().name[well_index]) {
return polymerFraction;
}
}
}
// we didn't find it return 0;
assert(false);
return 0.0;
}
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const
{
if (global_cell) {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
}
else {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(i, i));
}
}
}
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
computeRepRadiusPerfLength(const auto& grid)
{
// TODO, the function does not work for parallel running
// to be fixed later.
int number_of_cells = Opm::UgGridHelpers::numCells(grid);
const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
if (wells_ecl_.size() == 0) {
OPM_MESSAGE("No wells specified in Schedule section, "
"initializing no wells");
return;
}
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const size_t timeStep = current_timeIdx_;
wells_rep_radius_.clear();
wells_perf_length_.clear();
wells_bore_diameter_.clear();
wells_rep_radius_.reserve(nperf);
wells_perf_length_.reserve(nperf);
wells_bore_diameter_.reserve(nperf);
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(global_cell, number_of_cells,
cartesian_to_compressed);
int well_index = 0;
for (auto wellIter= wells_ecl_.begin(); wellIter != wells_ecl_.end(); ++wellIter) {
const auto* well = (*wellIter);
if (well->getStatus(timeStep) == WellCommon::SHUT) {
continue;
}
{ // COMPDAT handling
const auto& completionSet = well->getCompletions(timeStep);
for (size_t c=0; c<completionSet.size(); c++) {
const auto& completion = completionSet.get(c);
if (completion.getState() == WellCompletion::OPEN) {
int i = completion.getI();
int j = completion.getJ();
int k = completion.getK();
const int* cpgdim = cart_dims;
int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << well->name() << ')');
}
int cell = cgit->second;
{
double radius = 0.5*completion.getDiameter();
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
WellCompletion::DirectionEnum direction = completion.getDirection();
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
break;
case Opm::WellCompletion::DirectionEnum::Y:
re = std::sqrt(cubical[0] * cubical[2] / M_PI);
perf_length = cubical[1];
break;
case Opm::WellCompletion::DirectionEnum::Z:
re = std::sqrt(cubical[0] * cubical[1] / M_PI);
perf_length = cubical[2];
break;
default:
OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
}
double repR = std::sqrt(re * radius);
wells_rep_radius_.push_back(repR);
wells_perf_length_.push_back(perf_length);
wells_bore_diameter_.push_back(2. * radius);
}
} else {
if (completion.getState() != WellCompletion::SHUT) {
OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
}
}
}
}
well_index++;
}
}
} // namespace Opm