OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 10:40:21 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1279 from GitPaean/multisegment_well_refactoring_update

Browse Source 
Support for multi-segment wells
This commit is contained in:
Atgeirr Flø Rasmussen 2017-10-16 22:46:21 +02:00 committed by GitHub
commit 1858b56b07
18 changed files with 3060 additions and 322 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists_files.cmake
View File

@ -250,7 +250,11 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/WellInterface_impl.hpp
opm/autodiff/StandardWell.hpp
opm/autodiff/StandardWell_impl.hpp
opm/autodiff/StandardWellsDense.hpp
opm/autodiff/MultisegmentWell.hpp
opm/autodiff/MultisegmentWell_impl.hpp
opm/autodiff/MSWellHelpers.hpp
opm/autodiff/BlackoilWellModel.hpp
opm/autodiff/BlackoilWellModel_impl.hpp
opm/autodiff/StandardWellsSolvent.hpp
opm/autodiff/StandardWellsSolvent_impl.hpp
opm/autodiff/MissingFeatures.hpp

37
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -28,7 +28,7 @@
#include <ewoms/common/start.hh>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
@ -154,7 +154,7 @@ namespace Opm {
/// \param[in] terminal_output request output to cout/cerr
BlackoilModelEbos(Simulator& ebosSimulator,
const ModelParameters& param,
StandardWellsDense<TypeTag>& well_model,
BlackoilWellModel<TypeTag>& well_model,
RateConverterType& rate_converter,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output
@ -391,7 +391,7 @@ namespace Opm {
}
catch ( const Dune::FMatrixError& e )
{
OPM_THROW(Opm::NumericalProblem,"Well equation did not converge");
OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations");
}
return report;
@ -425,12 +425,12 @@ namespace Opm {
saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
}
@ -446,16 +446,16 @@ namespace Opm {
saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
}
Scalar tmp = pressureNew - pressureOld;
resultDelta += tmp*tmp;
resultDenom += pressureNew*pressureNew;
@ -503,14 +503,14 @@ namespace Opm {
// Solve system.
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense<TypeTag>, true > Operator;
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
Operator opA(ebosJac, well_model_, istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, StandardWellsDense<TypeTag>, false > Operator;
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, well_model_);
istlSolver().solve( opA, x, ebosResid );
}
@ -680,7 +680,7 @@ namespace Opm {
maxVal = std::max(std::abs(dsw),maxVal);
maxVal = std::max(std::abs(dsg),maxVal);
maxVal = std::max(std::abs(dso),maxVal);
maxVal = std::max(std::abs(dss),maxVal);
maxVal = std::max(std::abs(dss),maxVal);
double satScaleFactor = 1.0;
if (maxVal > dsMax()) {
@ -905,7 +905,6 @@ namespace Opm {
bool converged_MB = true;
bool converged_CNV = true;
bool converged_Well = true;
// Finish computation
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
@ -913,13 +912,12 @@ namespace Opm {
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
// Well flux convergence is only for fluid phases, not other materials
// in our current implementation.
converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
residual_norms.push_back(CNV[compIdx]);
}
const bool converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
bool converged = converged_MB && converged_Well;
// do not care about the cell based residual in the last two Newton
@ -1502,7 +1500,7 @@ namespace Opm {
SimulatorReport failureReport_;
// Well Model
StandardWellsDense<TypeTag>& well_model_;
BlackoilWellModel<TypeTag>& well_model_;
/// \brief Whether we print something to std::cout
bool terminal_output_;
@ -1519,9 +1517,10 @@ namespace Opm {
public:
/// return the StandardWells object
StandardWellsDense<TypeTag>&
BlackoilWellModel<TypeTag>&
wellModel() { return well_model_; }
const StandardWellsDense<TypeTag>&
const BlackoilWellModel<TypeTag>&
wellModel() const { return well_model_; }
int numWells() const { return well_model_.numWells(); }

14
opm/autodiff/BlackoilModelParameters.cpp
View File

@ -50,12 +50,20 @@ namespace Opm
tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_);
tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
tolerance_well_control_ = param.getDefault("tolerance_well_control", tolerance_well_control_);
max_welleq_iter_ = param.getDefault("max_welleq_iter", max_welleq_iter_);
if (use_multisegment_well_) {
tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_);
max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_);
use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_);
max_inner_iter_ms_wells_ = param.getDefault("max_inner_iter_ms_wells", max_inner_iter_ms_wells_);
}
maxSinglePrecisionTimeStep_ = unit::convert::from(
param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day );
max_strict_iter_ = param.getDefault("max_strict_iter",8);
solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
use_multisegment_well_ = param.getDefault("use_multisegment_well", use_multisegment_well_);
deck_file_name_ = param.template get<std::string>("deck_filename");
}
@ -75,10 +83,16 @@ namespace Opm
tolerance_cnv_ = 1.0e-2;
tolerance_wells_ = 1.0e-4;
tolerance_well_control_ = 1.0e-7;
tolerance_pressure_ms_wells_ = unit::convert::from(0.01, unit::barsa); // 0.01 bar
max_welleq_iter_ = 15;
max_pressure_change_ms_wells_ = unit::convert::from(2.0, unit::barsa); // 2.0 bar
use_inner_iterations_ms_wells_ = true;
max_inner_iter_ms_wells_ = 10;
maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day );
solve_welleq_initially_ = true;
update_equations_scaling_ = false;
use_update_stabilization_ = true;
use_multisegment_well_ = false;
}

22
opm/autodiff/BlackoilModelParameters.hpp
View File

@ -51,6 +51,19 @@ namespace Opm
/// Tolerance for the well control equations
// TODO: it might need to distinguish between rate control and pressure control later
double tolerance_well_control_;
/// Tolerance for the pressure equations for multisegment wells
double tolerance_pressure_ms_wells_;
/// Maximum pressure change over an iteratio for ms wells
double max_pressure_change_ms_wells_;
/// Whether to use inner iterations for ms wells
bool use_inner_iterations_ms_wells_;
/// Maximum inner iteration number for ms wells
int max_inner_iter_ms_wells_;
/// Maximum iteration number of the well equation solution
int max_welleq_iter_;
/// Tolerance for time step in seconds where single precision can be used
/// for solving for the Jacobian
@ -68,7 +81,14 @@ namespace Opm
/// Try to detect oscillation or stagnation.
bool use_update_stabilization_;
// The file name of the deck
/// Whether to use MultisegmentWell to handle multisegment wells
/// it is something temporary before the multisegment well model is considered to be
/// well developed and tested.
/// if it is false, we will handle multisegment wells as standard wells, which will be
/// the default behavoir for the moment. Later, we might set it to be true by default if necessary
bool use_multisegment_well_;
/// The file name of the deck
std::string deck_file_name_;
/// Construct from user parameters or defaults.

57
opm/autodiff/StandardWellsDense.hpp → opm/autodiff/BlackoilWellModel.hpp
View File

@ -21,8 +21,8 @@
*/
#ifndef OPM_STANDARDWELLSDENSE_HEADER_INCLUDED
#define OPM_STANDARDWELLSDENSE_HEADER_INCLUDED
#ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
#define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
#include <opm/common/OpmLog/OpmLog.hpp>
@ -49,6 +49,7 @@
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/WellInterface.hpp>
#include <opm/autodiff/StandardWell.hpp>
#include <opm/autodiff/MultisegmentWell.hpp>
#include<dune/common/fmatrix.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/matrixmatrix.hh>
@ -60,9 +61,9 @@
namespace Opm {
/// Class for handling the standard well model.
/// Class for handling the blackoil well model.
template<typename TypeTag>
class StandardWellsDense {
class BlackoilWellModel {
public:
// --------- Types ---------
typedef WellStateFullyImplicitBlackoil WellState;
@ -90,14 +91,14 @@ namespace Opm {
SurfaceToReservoirVoidage<BlackoilPropsAdFromDeck::FluidSystem, std::vector<int> >;
// --------- Public methods ---------
StandardWellsDense(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const RateConverterType& rate_converter,
const bool terminal_output,
const int current_index,
std::vector<int>& pvt_region_idx);
BlackoilWellModel(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const RateConverterType& rate_converter,
const bool terminal_output,
const int current_index,
const std::vector<int>& pvt_region_idx);
void init(const PhaseUsage phase_usage_arg,
const std::vector<bool>& active_arg,
@ -137,7 +138,7 @@ namespace Opm {
/// return true if wells are available on this process
bool localWellsActive() const;
bool getWellConvergence(Simulator& ebosSimulator,
bool getWellConvergence(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const;
/// upate the dynamic lists related to economic limits
@ -162,6 +163,8 @@ namespace Opm {
const int number_of_phases_;
const ModelParameters& param_;
using WellInterfacePtr = std::unique_ptr<WellInterface<TypeTag> >;
// a vector of all the wells.
// eventually, the wells_ above should be gone.
@ -174,12 +177,13 @@ namespace Opm {
// create the well container
static std::vector<WellInterfacePtr > createWellContainer(const Wells* wells,
const std::vector<const Well*>& wells_ecl,
const int time_step);
const bool use_multisegment_well,
const int time_step,
const ModelParameters& param);
// Well collection is used to enforce the group control
WellCollection* well_collection_;
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;
@ -188,7 +192,7 @@ namespace Opm {
PhaseUsage phase_usage_;
std::vector<bool> active_;
const RateConverterType& rate_converter_;
std::vector<int> pvt_region_idx_;
const std::vector<int>& pvt_region_idx_;
// the number of the cells in the local grid
int number_of_cells_;
@ -210,7 +214,7 @@ namespace Opm {
void computeRepRadiusPerfLength(const Grid& grid);
void computeAverageFormationFactor(Simulator& ebosSimulator,
void computeAverageFormationFactor(const Simulator& ebosSimulator,
std::vector<double>& B_avg) const;
void applyVREPGroupControl(WellState& well_state) const;
@ -229,15 +233,18 @@ namespace Opm {
void calculateEfficiencyFactors();
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw) const;
// it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation()
// makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions
// twice at the beginning of the time step
/// Calculating the explict quantities used in the well calculation. By explicit, we mean they are cacluated
/// at the beginning of the time step and no derivatives are included in these quantities
void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& xw) const;
SimulatorReport solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state) const;
void computeAccumWells() const;
void initPrimaryVariablesEvaluation() const;
// The number of components in the model.
@ -268,11 +275,15 @@ namespace Opm {
// some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
void prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state);
WellState& well_state) const;
void prepareGroupControl(const Simulator& ebos_simulator,
WellState& well_state) const;
};
} // namespace Opm
#include "StandardWellsDense_impl.hpp"
#include "BlackoilWellModel_impl.hpp"
#endif

226
opm/autodiff/StandardWellsDense_impl.hpp → opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -4,23 +4,23 @@
namespace Opm {
template<typename TypeTag>
StandardWellsDense<TypeTag>::
StandardWellsDense(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const RateConverterType& rate_converter,
const bool terminal_output,
const int current_timeIdx,
std::vector<int>& pvt_region_idx)
BlackoilWellModel<TypeTag>::
BlackoilWellModel(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const ModelParameters& param,
const RateConverterType& rate_converter,
const bool terminal_output,
const int current_timeIdx,
const std::vector<int>& pvt_region_idx)
: wells_active_(wells_arg!=nullptr)
, wells_(wells_arg)
, wells_ecl_(wells_ecl)
, number_of_wells_(wells_arg ? (wells_arg->number_of_wells) : 0)
, number_of_phases_(wells_arg ? (wells_arg->number_of_phases) : 0) // TODO: not sure if it is proper for this way
, well_container_(createWellContainer(wells_arg, wells_ecl, current_timeIdx) )
, well_collection_(well_collection)
, param_(param)
, well_container_(createWellContainer(wells_arg, wells_ecl, param.use_multisegment_well_, current_timeIdx, param) )
, well_collection_(well_collection)
, terminal_output_(terminal_output)
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
@ -36,7 +36,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
init(const PhaseUsage phase_usage_arg,
const std::vector<bool>& active_arg,
const double gravity_arg,
@ -87,7 +87,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
setVFPProperties(const VFPProperties* vfp_properties_arg)
{
for (auto& well : well_container_) {
@ -102,7 +102,7 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
numWells() const
{
return number_of_wells_;
@ -113,11 +113,13 @@ namespace Opm {
template<typename TypeTag>
std::vector<typename StandardWellsDense<TypeTag>::WellInterfacePtr >
StandardWellsDense<TypeTag>::
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>::
createWellContainer(const Wells* wells,
const std::vector< const Well* >& wells_ecl,
const int time_step)
const bool use_multisegment_well,
const int time_step,
const ModelParameters& param)
{
std::vector<WellInterfacePtr> well_container;
@ -146,13 +148,12 @@ namespace Opm {
}
const Well* well_ecl = wells_ecl[index_well];
// TODO: stopping throwing when encoutnering MS wells for now.
/* if (well_ecl->isMultiSegment(time_step)) {
OPM_THROW(Opm::NumericalProblem, "Not handling Multisegment Wells for now");
} */
// Basically, we are handling all the wells as StandardWell for the moment
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells) );
if ( !well_ecl->isMultiSegment(time_step) || !use_multisegment_well) {
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells, param) );
} else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells, param) );
}
}
}
return well_container;
@ -164,29 +165,27 @@ namespace Opm {
template<typename TypeTag>
SimulatorReport
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state)
{
if (iterationIdx == 0) {
prepareTimeStep(ebosSimulator, well_state);
}
SimulatorReport report;
if ( ! wellsActive() ) {
return report;
}
if (iterationIdx == 0) {
prepareTimeStep(ebosSimulator, well_state);
}
updateWellControls(well_state);
// Set the well primary variables based on the value of well solutions
initPrimaryVariablesEvaluation();
if (iterationIdx == 0) {
computeWellConnectionPressures(ebosSimulator, well_state);
computeAccumWells();
calculateExplicitQuantities(ebosSimulator, well_state);
}
if (param_.solve_welleq_initially_ && iterationIdx == 0) {
@ -205,7 +204,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
@ -224,7 +223,7 @@ namespace Opm {
// r = r - duneC_^T * invDuneD_ * resWell_
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
apply( BVector& r) const
{
if ( ! localWellsActive() ) {
@ -243,7 +242,7 @@ namespace Opm {
// Ax = A x - C D^-1 B x
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
// TODO: do we still need localWellsActive()?
@ -263,7 +262,7 @@ namespace Opm {
// Ax = Ax - alpha * C D^-1 B x
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
{
if ( ! localWellsActive() ) {
@ -287,11 +286,11 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const
{
for (auto& well : well_container_) {
well->recoverWellSolutionAndUpdateWellState(x, param_, well_state);
well->recoverWellSolutionAndUpdateWellState(x, well_state);
}
}
@ -301,7 +300,7 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
const auto& pu = phase_usage_;
@ -323,7 +322,7 @@ namespace Opm {
template<typename TypeTag>
int
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
numPhases() const
{
return number_of_phases_;
@ -335,7 +334,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
resetWellControlFromState(const WellState& xw) const
{
const int nw = wells_->number_of_wells;
@ -351,7 +350,7 @@ namespace Opm {
template<typename TypeTag>
bool
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
wellsActive() const
{
return wells_active_;
@ -363,7 +362,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
setWellsActive(const bool wells_active)
{
wells_active_ = wells_active;
@ -375,7 +374,7 @@ namespace Opm {
template<typename TypeTag>
bool
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
localWellsActive() const
{
return number_of_wells_ > 0;
@ -387,7 +386,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
initPrimaryVariablesEvaluation() const
{
for (auto& well : well_container_) {
@ -399,23 +398,9 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
computeAccumWells() const
{
for (auto& well : well_container_) {
well->computeAccumWell();
}
}
template<typename TypeTag>
SimulatorReport
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state) const
@ -427,6 +412,8 @@ namespace Opm {
std::vector< Scalar > B_avg( numComp, Scalar() );
computeAverageFormationFactor(ebosSimulator, B_avg);
const int max_iter = param_.max_welleq_iter_;
int it = 0;
bool converged;
do {
@ -447,7 +434,7 @@ namespace Opm {
if( localWellsActive() )
{
for (auto& well : well_container_) {
well->solveEqAndUpdateWellState(param_, well_state);
well->solveEqAndUpdateWellState(well_state);
}
}
// updateWellControls uses communication
@ -458,7 +445,7 @@ namespace Opm {
updateWellControls(well_state);
initPrimaryVariablesEvaluation();
}
} while (it < 15);
} while (it < max_iter);
if (converged) {
if ( terminal_output_ ) {
@ -490,14 +477,14 @@ namespace Opm {
template<typename TypeTag>
bool
StandardWellsDense<TypeTag>::
getWellConvergence(Simulator& ebosSimulator,
BlackoilWellModel<TypeTag>::
getWellConvergence(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const
{
ConvergenceReport report;
for (const auto& well : well_container_) {
report += well->getWellConvergence(ebosSimulator, B_avg, param_);
report += well->getWellConvergence(B_avg);
}
// checking NaN residuals
@ -549,14 +536,12 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw) const
BlackoilWellModel<TypeTag>::
calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& xw) const
{
if( ! localWellsActive() ) return ;
for (auto& well : well_container_) {
well->computeWellConnectionPressures(ebosSimulator, xw);
well->calculateExplicitQuantities(ebosSimulator, xw);
}
}
@ -566,7 +551,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
updateWellControls(WellState& xw) const
{
// Even if there no wells active locally, we cannot
@ -590,7 +575,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
updateListEconLimited(const Schedule& schedule,
const int current_step,
const Wells* wells_struct,
@ -608,18 +593,11 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const
{
updatePrimaryVariables(well_state);
computeWellConnectionPressures(ebosSimulator, well_state);
// initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
// TODO: for computeWellPotentials, no derivative is required actually
initPrimaryVariablesEvaluation();
// number of wells and phases
const int nw = number_of_wells_;
const int np = number_of_phases_;
@ -642,18 +620,50 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state)
WellState& well_state) const
{
const int nw = number_of_wells_;
for (int w = 0; w < nw; ++w) {
// after restarting, the well_controls can be modified while
// the well_state still uses the old control index
// we need to synchronize these two.
resetWellControlFromState(well_state);
// after restarting, the well_controls can be modified while
// the well_state still uses the old control index
// we need to synchronize these two.
// keep in mind that we set the control index of well_state to be the same with
// with the wellControls from the deck when we create well_state at the beginning of the report step
resetWellControlFromState(well_state);
if (wellCollection()->groupControlActive()) {
// process group control related
prepareGroupControl(ebos_simulator, well_state);
// since the controls are all updated, we should update well_state accordingly
for (int w = 0; w < number_of_wells_; ++w) {
WellControls* wc = well_container_[w]->wellControls();
const int control = well_controls_get_current(wc);
well_state.currentControls()[w] = control;
// TODO: for VFP control, the perf_densities are still zero here, investigate better
// way to handle it later.
well_container_[w]->updateWellStateWithTarget(control, well_state);
// The wells are not considered to be newly added
// for next time step
if (well_state.isNewWell(w) ) {
well_state.setNewWell(w, false);
}
} // end of for (int w = 0; w < nw; ++w)
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
prepareGroupControl(const Simulator& ebos_simulator,
WellState& well_state) const
{
// group control related processing
if (well_collection_->groupControlActive()) {
for (int w = 0; w < number_of_wells_; ++w) {
WellControls* wc = well_container_[w]->wellControls();
WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name());
@ -680,9 +690,7 @@ namespace Opm {
well_node.setIndividualControl(true);
}
}
}
if (well_collection_->groupControlActive()) {
if (well_collection_->requireWellPotentials()) {
// calculate the well potentials
@ -703,22 +711,6 @@ namespace Opm {
wellCollection()->updateWellTargets(well_state.wellRates());
}
}
// since the controls are all updated, we should update well_state accordingly
for (int w = 0; w < nw; ++w) {
WellControls* wc = well_container_[w]->wellControls();
const int control = well_controls_get_current(wc);
well_state.currentControls()[w] = control;
// TODO: for VFP control, the perf_densities are still zero here, investigate better
// way to handle it later.
well_container_[w]->updateWellStateWithTarget(control, well_state);
// The wells are not considered to be newly added
// for next time step
if (well_state.isNewWell(w) ) {
well_state.setNewWell(w, false);
}
} // end of for (int w = 0; w < nw; ++w)
}
@ -727,7 +719,7 @@ namespace Opm {
template<typename TypeTag>
WellCollection*
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
wellCollection() const
{
return well_collection_;
@ -739,7 +731,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
calculateEfficiencyFactors()
{
if ( !localWellsActive() ) {
@ -764,7 +756,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
computeWellVoidageRates(const WellState& well_state,
std::vector<double>& well_voidage_rates,
std::vector<double>& voidage_conversion_coeffs) const
@ -826,7 +818,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
applyVREPGroupControl(WellState& well_state) const
{
if ( wellCollection()->havingVREPGroups() ) {
@ -854,7 +846,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
updateGroupControls(WellState& well_state) const
{
@ -897,7 +889,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const
{
if (global_cell) {
@ -915,7 +907,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
computeRepRadiusPerfLength(const Grid& grid)
{
// TODO, the function does not work for parallel running
@ -937,8 +929,8 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
computeAverageFormationFactor(Simulator& ebosSimulator,
BlackoilWellModel<TypeTag>::
computeAverageFormationFactor(const Simulator& ebosSimulator,
std::vector<double>& B_avg) const
{
const int np = numPhases();
@ -984,7 +976,7 @@ namespace Opm {
template<typename TypeTag>
void
StandardWellsDense<TypeTag>::
BlackoilWellModel<TypeTag>::
updatePrimaryVariables(const WellState& well_state) const
{
for (const auto& well : well_container_) {

7
opm/autodiff/Compat.hpp
View File

@ -31,7 +31,6 @@ namespace Opm {
// Forward declarations
class SimulationDataContainer;
class WellStateFullyImplicitBlackoil;
class WellStateFullyImplicitBlackoilDense;
/// Extract single data vector from striped data.
/// \return u such that u[i] = v[offset + i*stride].
@ -71,12 +70,6 @@ namespace Opm {
PhaseUsage phases,
WellStateFullyImplicitBlackoil& state );
/// As the WellStateFullyImplicitBlackoil overload, but also sets
/// the wellSolution field from the values of the other fields.
void wellsToState( const data::Wells& wells,
PhaseUsage phases,
WellStateFullyImplicitBlackoilDense& state );
}
#endif //OPM_SIMULATORS_COMPAT_HPP

195
opm/autodiff/MSWellHelpers.hpp Normal file
View File

@ -0,0 +1,195 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
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_MSWELLHELPERS_HEADER_INCLUDED
#define OPM_MSWELLHELPERS_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>
#include <cmath>
namespace Opm {
namespace mswellhelpers
{
// obtain y = D^-1 * x with a direct solver
template <typename MatrixType, typename VectorType>
VectorType
invDXDirect(const MatrixType& D, VectorType x)
{
VectorType y(x.size());
y = 0.;
Dune::UMFPack<MatrixType> linsolver(D, 0);
// Object storing some statistics about the solving process
Dune::InverseOperatorResult res;
// Solve
linsolver.apply(y, x, res);
// Checking if there is any inf or nan in y
// it will be the solution before we find a way to catch the singularity of the matrix
for (size_t i_block = 0; i_block < y.size(); ++i_block) {
for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) {
if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) {
OPM_THROW(Opm::NumericalProblem, "nan or inf value found in invDXDirect due to singular matrix");
}
}
}
return y;
}
// obtain y = D^-1 * x with a BICSSTAB iterative solver
template <typename MatrixType, typename VectorType>
VectorType
invDX(const MatrixType& D, VectorType x)
{
// the function will change the value of x, so we should not use reference of x here.
// TODO: store some of the following information to avoid to call it again and again for
// efficiency improvement.
// Bassically, only the solve / apply step is different.
VectorType y(x.size());
y = 0.;
Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
// Sequential incomplete LU decomposition as the preconditioner
Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
// Dune::SeqILUn<MatrixType, VectorType, VectorType> preconditioner(D, 1, 0.92);
// Dune::SeqGS<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
// Dune::SeqJac<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
// Preconditioned BICGSTAB solver
Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
preconditioner,
1.e-8, // desired residual reduction factor
250, // maximum number of iterations
0); // verbosity of the solver */
// Object storing some statistics about the solving process
Dune::InverseOperatorResult res;
// Solve
linsolver.apply(y, x, res);
if ( !res.converged ) {
OPM_THROW(Opm::NumericalProblem, "the invDX does not get converged! ");
}
return y;
}
inline double haalandFormular(const double re, const double diameter, const double roughness)
{
const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
// sqrt(1/f) should be non-positive
assert(value >= 0.0);
return 1. / (value * value);
}
inline double calculateFrictionFactor(const double area, const double diameter,
const double w, const double roughness, const double mu)
{
double f = 0.;
// Reynolds number
const double re = std::abs(diameter * w / (area * mu));
if ( re == 0.0 ) {
// make sure it is because the mass rate is zero
assert(w == 0.);
return 0.0;
}
const double re_value1 = 200.;
const double re_value2 = 4000.;
if (re < re_value1) {
f = 16. / re;
} else if (re > re_value2){
f = haalandFormular(re, diameter, roughness);
} else { // in between
const double f1 = 16. / re_value1;
const double f2 = haalandFormular(re_value2, diameter, roughness);
f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
}
return f;
}
// calculating the friction pressure loss
// l is the segment length
// area is the segment cross area
// diameter is the segment inner diameter
// w is mass flow rate through the segment
// density is density
// roughness is the absolute roughness
// mu is the average phase viscosity
template <typename ValueType>
ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness,
const ValueType& density, const ValueType& w, const ValueType& mu)
{
const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value());
// \Note: a factor of 2 needs to be here based on the dimensional analysis
return 2. * f * l * w * w / (area * area * diameter * density);
}
template <typename ValueType>
ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density)
{
return (0.5 * mass_rate * mass_rate / (area * area * density));
}
} // namespace mswellhelpers
}
#endif

350
opm/autodiff/MultisegmentWell.hpp Normal file
View File

@ -0,0 +1,350 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
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_MULTISEGMENTWELL_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_HEADER_INCLUDED
#include <opm/autodiff/WellInterface.hpp>
namespace Opm
{
template<typename TypeTag>
class MultisegmentWell: public WellInterface<TypeTag>
{
public:
typedef WellInterface<TypeTag> Base;
using typename Base::WellState;
using typename Base::Simulator;
using typename Base::IntensiveQuantities;
using typename Base::FluidSystem;
using typename Base::ModelParameters;
using typename Base::MaterialLaw;
using typename Base::BlackoilIndices;
/// the number of reservior equations
using Base::numEq;
using Base::has_solvent;
using Base::has_polymer;
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
// TODO: the following system looks not rather flexible. Looking into all kinds of possibilities
// TODO: gas is always there? how about oil water case?
// Is it gas oil two phase case?
static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0);
static const int GTotal = 0;
static const int WFrac = gasoil? -1000: 1;
static const int GFrac = gasoil? 1 : 2;
static const int SPres = gasoil? 2 : 3;
/// the number of well equations // TODO: it should have a more general strategy for it
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq : numEq + 1;
using typename Base::Scalar;
using typename Base::ConvergenceReport;
/// the matrix and vector types for the reservoir
using typename Base::Mat;
using typename Base::BVector;
using typename Base::Eval;
// sparsity pattern for the matrices
// [A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
typedef Dune::FieldVector<Scalar, numWellEq> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// the matrix type for the diagonal matrix D
typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune::FieldMatrix<Scalar, numWellEq, numEq> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
// TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
// EvalR (Eval), EvalW, EvalRW
// TODO: for now, we only use one type to save some implementation efforts, while improve later.
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells);
virtual void initPrimaryVariablesEvaluation() const;
virtual void assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells);
/// updating the well state based the control mode specified with current
// TODO: later will check wheter we need current
virtual void updateWellStateWithTarget(const int current,
WellState& well_state) const;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const;
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const;
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const;
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials);
virtual void updatePrimaryVariables(const WellState& well_state) const;
virtual void solveEqAndUpdateWellState(WellState& well_state); // const?
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
/// number of segments for this well
/// int number_of_segments_;
int numberOfSegments() const;
int numberOfPerforations() const;
protected:
int number_segments_;
// components of the pressure drop to be included
WellSegment::CompPressureDropEnum compPressureDrop() const;
// multi-phase flow model
WellSegment::MultiPhaseModelEnum multiphaseModel() const;
// get the SegmentSet from the well_ecl_
const SegmentSet& segmentSet() const;
// protected member variables from the Base class
using Base::well_ecl_;
using Base::number_of_perforations_; // TODO: can use well_ecl_?
using Base::current_step_;
using Base::index_of_well_;
using Base::number_of_phases_;
// TODO: the current implementation really relies on the order of the
// perforation does not change from the parser to Wells structure.
using Base::well_cells_;
using Base::param_;
using Base::well_index_;
using Base::well_type_;
using Base::first_perf_;
using Base::saturation_table_number_;
using Base::well_efficiency_factor_;
using Base::gravity_;
using Base::well_controls_;
using Base::perf_depth_;
// protected functions from the Base class
using Base::active;
using Base::phaseUsage;
using Base::name;
using Base::numComponents;
using Base::flowPhaseToEbosPhaseIdx;
using Base::flowPhaseToEbosCompIdx;
using Base::getAllowCrossFlow;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int> > segment_perforations_;
// the inlet segments for each segment. It is for convinience and efficiency reason
std::vector<std::vector<int> > segment_inlets_;
// segment number is an ID of the segment, it is specified in the deck
// get the loation of the segment with a segment number in the segmentSet
int segmentNumberToIndex(const int segment_number) const;
// TODO, the following should go to a class for computing purpose
// two off-diagonal matrices
mutable OffDiagMatWell duneB_;
mutable OffDiagMatWell duneC_;
// diagonal matrix for the well
mutable DiagMatWell duneD_;
// residuals of the well equations
mutable BVectorWell resWell_;
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
mutable std::vector<std::array<double, numWellEq> > primary_variables_;
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
mutable std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
// depth difference between perforations and the perforated grid cells
std::vector<double> cell_perforation_depth_diffs_;
// pressure correction due to the different depth of the perforation and
// center depth of the grid block
std::vector<double> cell_perforation_pressure_diffs_;
// depth difference between the segment and the peforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_segment_depth_diffs_;
// the intial component compistion of segments
std::vector<std::vector<double> > segment_comp_initial_;
// the densities of segment fluids
// we should not have this member variable
std::vector<EvalWell> segment_densities_;
// the viscosity of the segments
std::vector<EvalWell> segment_viscosities_;
// the mass rate of the segments
std::vector<EvalWell> segment_mass_rates_;
std::vector<double> segment_depth_diffs_;
void initMatrixAndVectors(const int num_cells) const;
// protected functions
// EvalWell getBhp(); this one should be something similar to getSegmentPressure();
// EvalWell getQs(); this one should be something similar to getSegmentRates()
// EValWell wellVolumeFractionScaled, wellVolumeFraction, wellSurfaceVolumeFraction ... these should have different names, and probably will be needed.
// bool crossFlowAllowed(const Simulator& ebosSimulator) const; probably will be needed
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
const bool inner_iteration,
WellState& well_state) const;
// initialize the segment rates with well rates
// when there is no more accurate way to initialize the segment rates, we initialize
// the segment rates based on well rates with a simple strategy
void initSegmentRatesWithWellRates(WellState& well_state) const;
// computing the accumulation term for later use in well mass equations
void computeInitialComposition();
// compute the pressure difference between the perforation and cell center
void computePerfCellPressDiffs(const Simulator& ebosSimulator);
// fraction value of the primary variables
// should we just use member variables to store them instead of calculating them again and again
EvalWell volumeFraction(const int seg, const int comp_idx) const;
// F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
EvalWell volumeFractionScaled(const int seg, const int comp_idx) const;
// basically Q_p / \sigma_p Q_p
EvalWell surfaceVolumeFraction(const int seg, const int comp_idx) const;
void computePerfRate(const IntensiveQuantities& int_quants,
const std::vector<EvalWell>& mob_perfcells,
const int seg,
const int perf,
const EvalWell& segment_pressure,
const bool& allow_cf,
std::vector<EvalWell>& cq_s) const;
// convert a Eval from reservoir to contain the derivative related to wells
EvalWell extendEval(const Eval& in) const;
// compute the fluid properties, such as densities, viscosities, and so on, in the segments
// They will be treated implicitly, so they need to be of Evaluation type
void computeSegmentFluidProperties(const Simulator& ebosSimulator);
EvalWell getSegmentPressure(const int seg) const;
EvalWell getSegmentRate(const int seg, const int comp_idx) const;
EvalWell getSegmentGTotal(const int seg) const;
// get the mobility for specific perforation
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob) const;
void assembleControlEq() const;
void assemblePressureEq(const int seg) const;
// hytrostatic pressure loss
EvalWell getHydroPressureLoss(const int seg) const;
// frictinal pressure loss
EvalWell getFrictionPressureLoss(const int seg) const;
void handleAccelerationPressureLoss(const int seg) const;
// handling the overshooting and undershooting of the fractions
void processFractions(const int seg) const;
void updateWellStateFromPrimaryVariables(WellState& well_state) const;
double scalingFactor(const int comp_idx) const;
bool frictionalPressureLossConsidered() const;
bool accelerationalPressureLossConsidered() const;
// TODO: try to make ebosSimulator const, as it should be
void iterateWellEquations(Simulator& ebosSimulator,
const double dt,
WellState& well_state);
void assembleWellEqWithoutIteration(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells);
};
}
#include "MultisegmentWell_impl.hpp"
#endif // OPM_MULTISEGMENTWELL_HEADER_INCLUDED

1893
opm/autodiff/MultisegmentWell_impl.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

22
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -27,7 +27,7 @@
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/SimFIBODetails.hpp>
#include <opm/autodiff/moduleVersion.hpp>
@ -66,7 +66,7 @@ public:
typedef BlackoilModelEbos<TypeTag> Model;
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolver<Model> Solver;
typedef StandardWellsDense<TypeTag> WellModel;
typedef BlackoilWellModel<TypeTag> WellModel;
typedef RateConverter::SurfaceToReservoirVoidage<FluidSystem, std::vector<int> > RateConverterType;
@ -301,6 +301,20 @@ public:
WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_,
rateConverter_, terminal_output_, timer.currentStepNum(), legacyCellPvtRegionIdx_);
// handling MS well related
if (model_param_.use_multisegment_well_) { // if we use MultisegmentWell model
for (const auto& well : wells_ecl) {
// TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well
// while keep the well shut. More accurately, we should check if the well exisits in the Wells
// structure here
if (well->isMultiSegment(timer.currentStepNum()) ) { // there is one well is MS well
well_state.initWellStateMSWell(wells, wells_ecl, timer.currentStepNum(), phaseUsage_, prev_well_state);
break;
}
}
}
auto solver = createSolver(well_model);
std::vector<std::vector<double>> currentFluidInPlace;
@ -766,9 +780,9 @@ protected:
return totals;
}
void outputTimestampFIP(SimulatorTimer& timer, const std::string version)
{
{
std::ostringstream ss;
boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d %b %Y");
ss.imbue(std::locale(std::locale::classic(), facet));

35
opm/autodiff/StandardWell.hpp
View File

@ -99,7 +99,7 @@ namespace Opm
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
StandardWell(const Well* well, const int time_step, const Wells* wells);
StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
@ -120,21 +120,8 @@ namespace Opm
virtual void updateWellStateWithTarget(const int current,
WellState& xw) const;
// TODO: this should go to the WellInterface, while updateWellStateWithTarget
// will need touch different types of well_state, we will see.
virtual void updateWellControl(WellState& xw,
wellhelpers::WellSwitchingLogger& logger) const;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const;
/// computing the accumulation term for later use in well mass equations
virtual void computeAccumWell();
virtual void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw);
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const;
@ -143,16 +130,20 @@ namespace Opm
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param,
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const;
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
std::vector<double>& well_potentials) /* const */;
virtual void updatePrimaryVariables(const WellState& well_state) const;
virtual void solveEqAndUpdateWellState(WellState& well_state);
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
protected:
// protected functions from the Base class
@ -170,8 +161,8 @@ namespace Opm
// protected member variables from the Base class
using Base::vfp_properties_;
using Base::gravity_;
using Base::param_;
using Base::well_efficiency_factor_;
using Base::phase_usage_;
using Base::first_perf_;
using Base::ref_depth_;
using Base::perf_depth_;
@ -240,7 +231,6 @@ namespace Opm
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const;
// calculate the properties for the well connections
@ -268,8 +258,11 @@ namespace Opm
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state);
// computing the accumulation term for later use in well mass equations
void computeAccumWell();
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw);
// TODO: to check whether all the paramters are required
void computePerfRate(const IntensiveQuantities& intQuants,

147
opm/autodiff/StandardWell_impl.hpp
View File

@ -24,8 +24,8 @@ namespace Opm
{
template<typename TypeTag>
StandardWell<TypeTag>::
StandardWell(const Well* well, const int time_step, const Wells* wells)
: Base(well, time_step, wells)
StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
: Base(well, time_step, wells, param)
, perf_densities_(number_of_perforations_)
, perf_pressure_diffs_(number_of_perforations_)
, primary_variables_(numWellEq, 0.0)
@ -73,7 +73,6 @@ namespace Opm
}
for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
// Add nonzeros for diagonal
for (int perf = 0 ; perf < number_of_perforations_; ++perf) {
const int cell_idx = well_cells_[perf];
row.insert(cell_idx);
@ -317,7 +316,7 @@ namespace Opm
StandardWell<TypeTag>::
wellVolumeFraction(const int compIdx) const
{
const auto pu = phaseUsage();
const auto& pu = phaseUsage();
if (active()[Water] && compIdx == pu.phase_pos[Water]) {
return primary_variables_evaluation_[WFrac];
}
@ -401,14 +400,14 @@ namespace Opm
for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
}
auto& fs = intQuants.fluidState();
const auto& fs = intQuants.fluidState();
EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
EvalWell rs = extendEval(fs.Rs());
EvalWell rv = extendEval(fs.Rv());
const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const EvalWell rs = extendEval(fs.Rs());
const EvalWell rv = extendEval(fs.Rv());
std::vector<EvalWell> b_perfcells_dense(numComp, 0.0);
for (int phase = 0; phase < np; ++phase) {
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
}
if (has_solvent) {
@ -416,8 +415,8 @@ namespace Opm
}
// Pressure drawdown (also used to determine direction of flow)
EvalWell well_pressure = bhp + cdp;
EvalWell drawdown = pressure - well_pressure;
const EvalWell well_pressure = bhp + cdp;
const EvalWell drawdown = pressure - well_pressure;
// producing perforations
if ( drawdown.value() > 0 ) {
@ -523,8 +522,10 @@ namespace Opm
const int np = number_of_phases_;
// clear all entries
duneB_ = 0.0;
duneC_ = 0.0;
if (!only_wells) {
duneB_ = 0.0;
duneC_ = 0.0;
}
invDuneD_ = 0.0;
resWell_ = 0.0;
@ -763,12 +764,11 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const
{
const int np = number_of_phases_;
const double dBHPLimit = param.dbhp_max_rel_;
const double dFLimit = param.dwell_fraction_max_;
const double dBHPLimit = param_.dbhp_max_rel_;
const double dFLimit = param_.dwell_fraction_max_;
const auto pu = phaseUsage();
const std::vector<double> xvar_well_old = primary_variables_;
@ -1114,69 +1114,6 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
updateWellControl(WellState& xw,
wellhelpers::WellSwitchingLogger& logger) const
{
const int np = number_of_phases_;
const int w = index_of_well_;
const int old_control_index = xw.currentControls()[w];
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
WellControls* wc = well_controls_;
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
// the current control index
int current = xw.currentControls()[w];
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (wellhelpers::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, well_type_, wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
well_controls_set_current( wc, current);
}
// the new well control indices after all the related updates,
const int updated_control_index = xw.currentControls()[w];
// checking whether control changed
if (updated_control_index != old_control_index) {
logger.wellSwitched(name(),
well_controls_iget_type(wc, old_control_index),
well_controls_iget_type(wc, updated_control_index));
}
if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) {
updateWellStateWithTarget(updated_control_index, xw);
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
@ -1190,7 +1127,7 @@ namespace Opm
const int nperf = number_of_perforations_;
// TODO: can make this a member?
const int numComp = numComponents();
const PhaseUsage& pu = *phase_usage_;
const PhaseUsage& pu = phaseUsage();
b_perf.resize(nperf*numComp);
surf_dens_perf.resize(nperf*numComp);
const int w = index_of_well_;
@ -1298,7 +1235,7 @@ namespace Opm
const int np = number_of_phases_;
const int nperf = number_of_perforations_;
const int num_comp = numComponents();
const PhaseUsage* phase_usage = phase_usage_;
const PhaseUsage& phase_usage = phaseUsage();
// 1. Compute the flow (in surface volume units for each
// component) exiting up the wellbore from each perforation,
@ -1326,8 +1263,8 @@ namespace Opm
// absolute values of the surface rates divided by their sum.
// Then compute volume ratios (formation factors) for each perforation.
// Finally compute densities for the segments associated with each perforation.
const int gaspos = phase_usage->phase_pos[BlackoilPhases::Vapour];
const int oilpos = phase_usage->phase_pos[BlackoilPhases::Liquid];
const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
std::vector<double> mix(num_comp,0.0);
std::vector<double> x(num_comp);
std::vector<double> surf_dens(num_comp);
@ -1426,9 +1363,7 @@ namespace Opm
template<typename TypeTag>
typename StandardWell<TypeTag>::ConvergenceReport
StandardWell<TypeTag>::
getWellConvergence(Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const
getWellConvergence(const std::vector<double>& B_avg) const
{
typedef double Scalar;
typedef std::vector< Scalar > Vector;
@ -1440,8 +1375,8 @@ namespace Opm
// For the polymer case, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == numComp) || has_polymer);
const double tol_wells = param.tolerance_wells_;
const double maxResidualAllowed = param.max_residual_allowed_;
const double tol_wells = param_.tolerance_wells_;
const double maxResidualAllowed = param_.max_residual_allowed_;
// TODO: it should be the number of numWellEq
// using numComp here for flow_ebos running 2p case.
@ -1553,15 +1488,28 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state)
solveEqAndUpdateWellState(WellState& well_state)
{
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
BVectorWell dx_well(1);
invDuneD_.mv(resWell_, dx_well);
updateWellState(dx_well, param, well_state);
updateWellState(dx_well, well_state);
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state)
{
computeWellConnectionPressures(ebosSimulator, well_state);
computeAccumWell();
}
@ -1573,6 +1521,8 @@ namespace Opm
StandardWell<TypeTag>::
computeAccumWell()
{
// TODO: it should be num_comp, while it also bring problem for
// the polymer case.
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
}
@ -1642,12 +1592,11 @@ namespace Opm
void
StandardWell<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x,
const ModelParameters& param,
WellState& well_state) const
{
BVectorWell xw(1);
recoverSolutionWell(x, xw);
updateWellState(xw, param, well_state);
updateWellState(xw, well_state);
}
@ -1723,7 +1672,7 @@ namespace Opm
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
const Opm::PhaseUsage& pu = *phase_usage_;
const Opm::PhaseUsage& pu = phaseUsage();
std::vector<double> rates(3, 0.0);
if (active()[ Water ]) {
@ -1794,10 +1743,16 @@ namespace Opm
StandardWell<TypeTag>::
computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const
std::vector<double>& well_potentials) // const
{
const int np = number_of_phases_;
updatePrimaryVariables(well_state);
computeWellConnectionPressures(ebosSimulator, well_state);
// initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
// TODO: for computeWellPotentials, no derivative is required actually
initPrimaryVariablesEvaluation();
const int np = number_of_phases_;
well_potentials.resize(np, 0.0);
// get the bhp value based on the bhp constraints

36
opm/autodiff/WellInterface.hpp
View File

@ -86,7 +86,7 @@ namespace Opm
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
/// Constructor
WellInterface(const Well* well, const int time_step, const Wells* wells);
WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
/// Virutal destructor
virtual ~WellInterface() {}
@ -144,12 +144,9 @@ namespace Opm
}
};
virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const = 0;
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const = 0;
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state) = 0;
virtual void solveEqAndUpdateWellState(WellState& well_state) = 0;
virtual void assembleWellEq(Simulator& ebosSimulator,
const double dt,
@ -165,7 +162,7 @@ namespace Opm
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param,
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const = 0;
/// Ax = Ax - C D^-1 B x
@ -174,25 +171,22 @@ namespace Opm
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const = 0;
// TODO: before we decide to put more information under mutable, this function is not const
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const = 0;
virtual void computeAccumWell() = 0;
// TODO: it should come with a different name
// for MS well, the definition is different and should not use this name anymore
virtual void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw) = 0;
std::vector<double>& well_potentials) = 0;
virtual void updateWellStateWithTarget(const int current,
WellState& xw) const = 0;
virtual void updateWellControl(WellState& xw,
wellhelpers::WellSwitchingLogger& logger) const = 0;
void updateWellControl(WellState& xw,
wellhelpers::WellSwitchingLogger& logger) const;
virtual void updatePrimaryVariables(const WellState& well_state) const = 0;
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& xw) = 0; // should be const?
protected:
// to indicate a invalid connection
@ -205,6 +199,9 @@ namespace Opm
// the index of well in Wells struct
int index_of_well_;
// simulation parameters
const ModelParameters& param_;
// well type
// INJECTOR or PRODUCER
enum WellType well_type_;
@ -217,21 +214,18 @@ namespace Opm
std::vector<double> comp_frac_;
// controls for this well
// TODO: later will check whehter to let it stay with pointer
struct WellControls* well_controls_;
// number of the perforations for this well
int number_of_perforations_;
// record the index of the first perforation
// TODO: it might not be needed if we refactor WellState to be a vector
// of states of individual well.
int first_perf_;
// well index for each perforation
std::vector<double> well_index_;
// TODO: it might should go to StandardWell
// depth for each perforation
std::vector<double> perf_depth_;
@ -273,7 +267,7 @@ namespace Opm
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
// TODO: it is dumplicated with StandardWellsDense
// TODO: it is dumplicated with BlackoilWellModel
int numComponents() const;
double wsolvent() const;

70
opm/autodiff/WellInterface_impl.hpp
View File

@ -25,9 +25,10 @@ namespace Opm
template<typename TypeTag>
WellInterface<TypeTag>::
WellInterface(const Well* well, const int time_step, const Wells* wells)
WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
: well_ecl_(well)
, current_step_(time_step)
, param_(param)
{
if (!well) {
OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface");
@ -401,13 +402,76 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::
updateWellControl(WellState& well_state,
wellhelpers::WellSwitchingLogger& logger) const
{
const int np = number_of_phases_;
const int w = index_of_well_;
const int old_control_index = well_state.currentControls()[w];
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
WellControls* wc = well_controls_;
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
// the current control index
int current = well_state.currentControls()[w];
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (wellhelpers::constraintBroken(
well_state.bhp(), well_state.thp(), well_state.wellRates(),
w, np, well_type_, wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
well_state.currentControls()[w] = ctrl_index;
current = well_state.currentControls()[w];
well_controls_set_current( wc, current);
}
// the new well control indices after all the related updates,
const int updated_control_index = well_state.currentControls()[w];
// checking whether control changed
if (updated_control_index != old_control_index) {
logger.wellSwitched(name(),
well_controls_iget_type(wc, old_control_index),
well_controls_iget_type(wc, updated_control_index));
}
if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) {
updateWellStateWithTarget(updated_control_index, well_state);
}
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
const Opm::PhaseUsage& pu = *phase_usage_;
const Opm::PhaseUsage& pu = phaseUsage();
const int np = number_of_phases_;
if (econ_production_limits.onMinOilRate()) {
@ -464,7 +528,7 @@ namespace Opm
double violation_extent = -1.0;
const int np = number_of_phases_;
const Opm::PhaseUsage& pu = *phase_usage_;
const Opm::PhaseUsage& pu = phaseUsage();
const int well_number = index_of_well_;
assert(active()[Oil]);

4
opm/autodiff/WellMultiSegment.cpp
View File

@ -60,7 +60,7 @@ namespace Opm
// we change the ID to location now for easier use later.
for (int i = 0; i < m_number_of_segments_; ++i) {
// The segment number for top segment is 0, the segment number of its outlet segment will be -1
m_outlet_segment_[i] = segment_set.numberToLocation(segment_set[i].outletSegment());
m_outlet_segment_[i] = segment_set.segmentNumberToIndex(segment_set[i].outletSegment());
m_segment_length_[i] = segment_set[i].totalLength();
m_segment_depth_[i] = segment_set[i].depth();
m_segment_internal_diameter_[i] = segment_set[i].internalDiameter();
@ -111,7 +111,7 @@ namespace Opm
int i_segment = completion_set.get(i).getSegmentNumber();
// using the location of the segment in the array as the segment number/id.
// TODO: it can be helpful for output or postprocessing if we can keep the original number.
i_segment = segment_set.numberToLocation(i_segment);
i_segment = segment_set.segmentNumberToIndex(i_segment);
m_segment_perforations_[i_segment].push_back(i);
temp_perf_depth[i] = completion_set.get(i).getCenterDepth();
}

248
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -26,6 +26,7 @@
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <vector>
#include <cassert>
@ -165,7 +166,7 @@ namespace Opm
}
}
}
// perfPressures
// perfPressures
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
@ -206,6 +207,17 @@ namespace Opm
}
}
}
{
// we need to create a trival segment related values to avoid there will be some
// multi-segment wells added later.
top_segment_index_.reserve(nw);
for (int w = 0; w < nw; ++w) {
top_segment_index_.push_back(w);
}
segpress_ = bhp();
segrates_ = wellRates();
}
}
template <class State>
@ -247,7 +259,6 @@ namespace Opm
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) );
}
}
@ -282,6 +293,197 @@ namespace Opm
}
/// init the MS well related.
template <typename PrevWellState>
void initWellStateMSWell(const Wells* wells, const std::vector<const Well*>& wells_ecl,
const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state)
{
// still using the order in wells
const int nw = wells->number_of_wells;
if (nw == 0) {
return;
}
top_segment_index_.clear();
top_segment_index_.reserve(nw);
segpress_.clear();
segpress_.reserve(nw);
segrates_.clear();
segrates_.reserve(nw * numPhases());
nseg_ = 0;
// in the init function, the well rates and perforation rates have been initialized or copied from prevState
// what we do here, is to set the segment rates and perforation rates
for (int w = 0; w < nw; ++w) {
const int nw_wells_ecl = wells_ecl.size();
int index_well_ecl = 0;
const std::string well_name(wells->name[w]);
for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
if (well_name == wells_ecl[index_well_ecl]->name()) {
break;
}
}
// It should be able to find in wells_ecl.
if (index_well_ecl == nw_wells_ecl) {
OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
}
const Well* well_ecl = wells_ecl[index_well_ecl];
top_segment_index_.push_back(nseg_);
if ( !well_ecl->isMultiSegment(time_step) ) { // not multi-segment well
nseg_ += 1;
segpress_.push_back(bhp()[w]);
const int np = numPhases();
for (int p = 0; p < np; ++p) {
segrates_.push_back(wellRates()[np * w + p]);
}
} else { // it is a multi-segment well
const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step);
// assuming the order of the perforations in well_ecl is the same with Wells
const CompletionSet& completion_set = well_ecl->getCompletions(time_step);
// number of segment for this single well
const int well_nseg = segment_set.numberSegment();
const int nperf = completion_set.size();
nseg_ += well_nseg;
// we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment
// that is why I think we should use a well model to initialize the WellState here
std::vector<std::vector<int>> segment_perforations(well_nseg);
for (int perf = 0; perf < nperf; ++perf) {
const Completion& completion = completion_set.get(perf);
const int segment_number = completion.getSegmentNumber();
const int segment_index = segment_set.segmentNumberToIndex(segment_number);
segment_perforations[segment_index].push_back(perf);
}
std::vector<std::vector<int>> segment_inlets(well_nseg);
for (int seg = 0; seg < well_nseg; ++seg) {
const Segment& segment = segment_set[seg];
const int segment_number = segment.segmentNumber();
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) {
const int segment_index = segment_set.segmentNumberToIndex(segment_number);
const int outlet_segment_index = segment_set.segmentNumberToIndex(outlet_segment_number);
segment_inlets[outlet_segment_index].push_back(segment_index);
}
}
// for the segrates_, now it becomes a recursive solution procedure.
{
const int np = numPhases();
const int start_perf = wells->well_connpos[w];
const int start_perf_next_well = wells->well_connpos[w + 1];
assert(nperf == (start_perf_next_well - start_perf)); // make sure the information from wells_ecl consistent with wells
if (pu.phase_used[Gas]) {
const int gaspos = pu.phase_pos[Gas];
// scale the phase rates for Gas to avoid too bad initial guess for gas fraction
// it will probably benefit the standard well too, while it needs to be justified
// TODO: to see if this strategy can benefit StandardWell too
// TODO: it might cause big problem for gas rate control or if there is a gas rate limit
// maybe the best way is to initialize the fractions first then get the rates
for (int perf = 0; perf < nperf; perf++) {
const int perf_pos = start_perf + perf;
perfPhaseRates()[np * perf_pos + gaspos] *= 100.;
}
}
const std::vector<double> perforation_rates(perfPhaseRates().begin() + np * start_perf,
perfPhaseRates().begin() + np * start_perf_next_well); // the perforation rates for this well
std::vector<double> segment_rates;
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, segment_rates);
std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
}
// for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment
// if there is no perforation associated with this segment, it uses the pressure from the outlet segment
// which requres the ordering is successful
// Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be
// improved during the solveWellEq process
{
// top segment is always the first one, and its pressure is the well bhp
segpress_.push_back(bhp()[w]);
const int top_segment = top_segment_index_[w];
const int start_perf = wells->well_connpos[w];
for (int seg = 1; seg < well_nseg; ++seg) {
if ( !segment_perforations[seg].empty() ) {
const int first_perf = segment_perforations[seg][0];
segpress_.push_back(perfPress()[start_perf + first_perf]);
} else {
// segpress_.push_back(bhp); // may not be a good decision
// using the outlet segment pressure // it needs the ordering is correct
const int outlet_seg = segment_set[seg].outletSegment();
segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
}
}
}
}
}
assert(int(segpress_.size()) == nseg_);
assert(int(segrates_.size()) == nseg_ * numPhases() );
if (!prev_well_state.wellMap().empty()) {
// copying MS well related
const auto& end = prev_well_state.wellMap().end();
const int np = numPhases();
for (int w = 0; w < nw; ++w) {
const std::string name( wells->name[w] );
const auto& it = prev_well_state.wellMap().find( name );
if (it != end) { // the well is found in the prev_well_state
// TODO: the well with same name can change a lot, like they might not have same number of segments
// we need to handle that later.
// for now, we just copy them.
const int old_index_well = (*it).second[0];
const int new_index_well = w;
const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well);
const int new_top_segmnet_index = topSegmentIndex(new_index_well);
int number_of_segment = 0;
// if it is the last well in list
if (new_index_well == int(top_segment_index_.size()) - 1) {
number_of_segment = nseg_ - new_top_segmnet_index;
} else {
number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index;
}
for (int i = 0; i < number_of_segment * np; ++i) {
segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i];
}
for (int i = 0; i < number_of_segment; ++i) {
segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i];
}
}
}
}
}
static void calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets, const std::vector<std::vector<int>>&segment_perforations,
const std::vector<double>& perforation_rates, const int np, const int segment, std::vector<double>& segment_rates)
{
// the rate of the segment equals to the sum of the contribution from the perforations and inlet segment rates.
// the first segment is always the top segment, its rates should be equal to the well rates.
assert(segment_inlets.size() == segment_perforations.size());
const int well_nseg = segment_inlets.size();
if (segment == 0) { // beginning the calculation
segment_rates.resize(np * well_nseg, 0.0);
}
// contributions from the perforations belong to this segment
for (const int& perf : segment_perforations[segment]) {
for (int p = 0; p < np; ++p) {
segment_rates[np * segment + p] += perforation_rates[np * perf + p];
}
}
for (const int& inlet_seg : segment_inlets[segment]) {
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates);
for (int p = 0; p < np; ++p) {
segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p];
}
}
}
bool isNewWell(const int w) const {
return is_new_well_[w];
}
@ -305,6 +507,38 @@ namespace Opm
return solvent_well_rate;
}
const std::vector<double>& segRates() const
{
return segrates_;
}
std::vector<double>& segRates()
{
return segrates_;
}
const std::vector<double>& segPress() const
{
return segpress_;
}
std::vector<double>& segPress()
{
return segpress_;
}
int numSegment() const
{
return nseg_;
}
int topSegmentIndex(const int w) const
{
assert(w < int(top_segment_index_.size()) );
return top_segment_index_[w];
}
private:
std::vector<double> perfphaserates_;
std::vector<int> current_controls_;
@ -315,6 +549,16 @@ namespace Opm
// will have very wrong compositions for production wells, will mostly cause
// problem with VFP interpolation
std::vector<bool> is_new_well_;
// MS well related
// for StandardWell, the number of segments will be one
std::vector<double> segrates_;
std::vector<double> segpress_;
// the index of the top segments, which is used to locate the
// multisegment well related information in WellState
std::vector<int> top_segment_index_;
int nseg_; // total number of the segments
};
} // namespace Opm

13
tests/test_wellmodel.cpp
View File

@ -47,6 +47,7 @@
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
@ -56,7 +57,7 @@
#include <ewoms/common/start.hh>
#include <opm/autodiff/StandardWell.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
// maybe should just include BlackoilModelEbos.hpp
namespace Ewoms {
@ -122,9 +123,10 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const auto& wells_ecl = setup_test.ecl_state->getSchedule().getWells(setup_test.current_timestep);
BOOST_CHECK_EQUAL( wells_ecl.size(), 2);
const Opm::Well* well = wells_ecl[1];
BOOST_CHECK_THROW( StandardWell( well, -1, wells), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr), std::invalid_argument);
const Opm::BlackoilModelParameters param;
BOOST_CHECK_THROW( StandardWell( well, -1, wells, param), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param), std::invalid_argument);
}
@ -137,6 +139,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{
const int nw = wells_struct ? (wells_struct->number_of_wells) : 0;
const Opm::BlackoilModelParameters param;
for (int w = 0; w < nw; ++w) {
const std::string well_name(wells_struct->name[w]);
@ -150,7 +153,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
// we should always be able to find the well in wells_ecl
BOOST_CHECK(index_well != wells_ecl.size());
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct) );
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param) );
}
}