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

removing StandardWellV

Browse Source
This commit is contained in:
Kai Bao 2019-03-25 12:52:34 +01:00
parent 33de77e595
commit 3dcec6a436
7 changed files with 421 additions and 3742 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -206,8 +206,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/StandardWell_impl.hpp
opm/simulators/wells/MultisegmentWell.hpp
opm/simulators/wells/MultisegmentWell_impl.hpp
opm/simulators/wells/StandardWellV.hpp
opm/simulators/wells/StandardWellV_impl.hpp
opm/simulators/wells/MSWellHelpers.hpp
opm/simulators/wells/BlackoilWellModel.hpp
opm/simulators/wells/BlackoilWellModel_impl.hpp

1
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -46,7 +46,6 @@
#include <opm/autodiff/RateConverter.hpp>
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/StandardWell.hpp>
#include <opm/simulators/wells/StandardWellV.hpp>
#include <opm/simulators/wells/MultisegmentWell.hpp>
#include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
#include <opm/autodiff/SimFIBODetails.hpp>

55
opm/simulators/wells/StandardWell.hpp
View File

@ -28,6 +28,11 @@
#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/material/densead/DynamicEvaluation.hpp>
#include <dune/common/dynvector.hh>
#include <dune/common/dynmatrix.hh>
namespace Opm
{
@ -89,10 +94,6 @@ namespace Opm
// TODO: we should have indices for the well equations and well primary variables separately
static const int Bhp = numStaticWellEq - numWellControlEq;
// total number of the well equations and primary variables
// for StandardWell, no extra well equations will be used.
static const int numWellEq = numStaticWellEq;
using typename Base::Scalar;
@ -110,18 +111,18 @@ namespace Opm
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
typedef Dune::FieldVector<Scalar, numWellEq> VectorBlockWellType;
typedef Dune::DynamicVector<Scalar> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// the matrix type for the diagonal matrix D l
typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
// the matrix type for the diagonal matrix D
typedef Dune::DynamicMatrix<Scalar> 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::DynamicMatrix<Scalar> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
typedef DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + numEq + 1> EvalWell;
using Base::contiSolventEqIdx;
using Base::contiPolymerEqIdx;
@ -229,6 +230,10 @@ namespace Opm
using Base::perf_length_;
using Base::bore_diameters_;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
// densities of the fluid in each perforation
std::vector<double> perf_densities_;
// pressure drop between different perforations
@ -380,15 +385,13 @@ namespace Opm
// mostly related to BHP limit and THP limit
virtual void checkWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
) override;
Opm::DeferredLogger& deferred_logger) override;
// check whether the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
void updateWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
);
Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under BHP limit with current reservoir condition
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
@ -444,6 +447,32 @@ namespace Opm
WellState& well_state, WellTestState& welltest_state,
Opm::DeferredLogger& deferred_logger) override;
// calculate the skin pressure based on water velocity, throughput and polymer concentration.
// throughput is used to describe the formation damage during water/polymer injection.
// calculated skin pressure will be applied to the drawdown during perforation rate calculation
// to handle the effect from formation damage.
EvalWell pskin(const double throuhgput,
const EvalWell& water_velocity,
const EvalWell& poly_inj_conc,
Opm::DeferredLogger& deferred_logger) const;
// calculate the skin pressure based on water velocity, throughput during water injection.
EvalWell pskinwater(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const;
// calculate the injecting polymer molecular weight based on the througput and water velocity
EvalWell wpolymermw(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const;
// handle the extra equations for polymer injectivity study
void handleInjectivityRateAndEquations(const IntensiveQuantities& int_quants,
const WellState& well_state,
const int perf,
std::vector<EvalWell>& cq_s,
Opm::DeferredLogger& deferred_logger);
virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
void checkConvergenceControlEq(ConvergenceReport& report,

483
opm/simulators/wells/StandardWellV.hpp
View File

@ -1,483 +0,0 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 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/>.
*/
#ifndef OPM_STANDARDWELLV_HEADER_INCLUDED
#define OPM_STANDARDWELLV_HEADER_INCLUDED
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/material/densead/DynamicEvaluation.hpp>
#include <dune/common/dynvector.hh>
#include <dune/common/dynmatrix.hh>
namespace Opm
{
template<typename TypeTag>
class StandardWellV: public WellInterface<TypeTag>
{
public:
typedef WellInterface<TypeTag> Base;
// TODO: some functions working with AD variables handles only with values (double) without
// dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value.
// And also, it can also be beneficial to make these functions hanle different types of AD variables.
using typename Base::Simulator;
using typename Base::WellState;
using typename Base::IntensiveQuantities;
using typename Base::FluidSystem;
using typename Base::MaterialLaw;
using typename Base::ModelParameters;
using typename Base::Indices;
using typename Base::PolymerModule;
using typename Base::RateConverterType;
using Base::numEq;
using Base::has_solvent;
using Base::has_polymer;
using Base::has_energy;
// polymer concentration and temperature are already known by the well, so
// polymer and energy conservation do not need to be considered explicitly
static const int numPolymerEq = Indices::numPolymers;
static const int numEnergyEq = Indices::numEnergy;
// number of the conservation equations
static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq;
// number of the well control equations
static const int numWellControlEq = 1;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
static const int numStaticWellEq = numWellConservationEq + numWellControlEq;
// the positions of the primary variables for StandardWell
// the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
// which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP.
// correspondingly, we have four well equations for blackoil model, the first three are mass
// converstation equations, and the last one is the well control equation.
// primary variables related to other components, will be before the Bhp and after F_g.
// well control equation is always the last well equation.
// TODO: in the current implementation, we use the well rate as the first primary variables for injectors,
// instead of G_t.
static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
static const int WQTotal = 0;
static const int WFrac = gasoil? -1000: 1;
static const int GFrac = gasoil? 1: 2;
static const int SFrac = !has_solvent ? -1000 : 3;
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
static const int Bhp = numStaticWellEq - numWellControlEq;
using typename Base::Scalar;
using Base::name;
using Base::Water;
using Base::Oil;
using Base::Gas;
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::DynamicVector<Scalar> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// the matrix type for the diagonal matrix D
typedef Dune::DynamicMatrix<Scalar> DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune::DynamicMatrix<Scalar> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
typedef DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + numEq + 1> EvalWell;
using Base::contiSolventEqIdx;
using Base::contiPolymerEqIdx;
static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
StandardWellV(const Well2& well, const int time_step, const Wells* wells,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells) override;
virtual void initPrimaryVariablesEvaluation() const override;
virtual void assembleWellEq(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double dt,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) override;
virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const override;
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const override;
/// 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,
Opm::DeferredLogger& deferred_logger) const override;
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) /* const */ override;
virtual void updatePrimaryVariables(const WellState& well_state, Opm::DeferredLogger& deferred_logger) const override;
virtual void solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger) override;
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override; // should be const?
virtual void addWellContributions(Mat& mat) const override;
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const override
{
return param_.matrix_add_well_contributions_;
}
protected:
// protected functions from the Base class
using Base::getAllowCrossFlow;
using Base::phaseUsage;
using Base::flowPhaseToEbosCompIdx;
using Base::ebosCompIdxToFlowCompIdx;
using Base::wsolvent;
using Base::wpolymer;
using Base::wellHasTHPConstraints;
using Base::mostStrictBhpFromBhpLimits;
using Base::scalingFactor;
using Base::scaleProductivityIndex;
// protected member variables from the Base class
using Base::current_step_;
using Base::well_ecl_;
using Base::vfp_properties_;
using Base::gravity_;
using Base::param_;
using Base::well_efficiency_factor_;
using Base::first_perf_;
using Base::ref_depth_;
using Base::perf_depth_;
using Base::well_cells_;
using Base::number_of_perforations_;
using Base::number_of_phases_;
using Base::saturation_table_number_;
using Base::comp_frac_;
using Base::well_index_;
using Base::index_of_well_;
using Base::well_controls_;
using Base::well_type_;
using Base::num_components_;
using Base::connectionRates_;
using Base::perf_rep_radius_;
using Base::perf_length_;
using Base::bore_diameters_;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
// densities of the fluid in each perforation
std::vector<double> perf_densities_;
// pressure drop between different perforations
std::vector<double> perf_pressure_diffs_;
// residuals of the well equations
BVectorWell resWell_;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
DiagMatWell invDuneD_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
mutable std::vector<double> primary_variables_;
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
mutable std::vector<EvalWell> primary_variables_evaluation_;
// the saturations in the well bore under surface conditions at the beginning of the time step
std::vector<double> F0_;
// the vectors used to describe the inflow performance relationship (IPR)
// Q = IPR_A - BHP * IPR_B
// TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
// it is only updated and used for producers for now
mutable std::vector<double> ipr_a_;
mutable std::vector<double> ipr_b_;
const EvalWell& getBhp() const;
EvalWell getQs(const int comp_idx) const;
const EvalWell& getWQTotal() const;
EvalWell wellVolumeFractionScaled(const int phase) const;
EvalWell wellVolumeFraction(const unsigned compIdx) const;
EvalWell wellSurfaceVolumeFraction(const int phase) const;
EvalWell extendEval(const Eval& in) const;
// 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,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
// calculate the properties for the well connections
// to calulate the pressure difference between well connections.
void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf) const;
// TODO: not total sure whether it is a good idea to put this function here
// the major reason to put here is to avoid the usage of Wells struct
void computeConnectionDensities(const std::vector<double>& perfComponentRates,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
void computeConnectionPressureDelta();
void computeWellConnectionDensitesPressures(const WellState& well_state,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
// computing the accumulation term for later use in well mass equations
void computeAccumWell();
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state);
void computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob,
const EvalWell& bhp,
const int perf,
const bool allow_cf,
std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger);
std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger);
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const;
double calculateThpFromBhp(const std::vector<double>& rates, const double bhp, Opm::DeferredLogger& deferred_logger) const;
// get the mobility for specific perforation
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob,
Opm::DeferredLogger& deferred_logger) const;
void updateWaterMobilityWithPolymer(const Simulator& ebos_simulator,
const int perf,
std::vector<EvalWell>& mob_water,
Opm::DeferredLogger& deferred_logger) const;
void updatePrimaryVariablesNewton(const BVectorWell& dwells,
const WellState& well_state) const;
void updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const;
void updateThp(WellState& well_state, Opm::DeferredLogger& deferred_logger) const;
void assembleControlEq(Opm::DeferredLogger& deferred_logger);
// handle the non reasonable fractions due to numerical overshoot
void processFractions() const;
// updating the inflow based on the current reservoir condition
void updateIPR(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) const;
// update the operability status of the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
virtual void checkWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override;
// check whether the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
void updateWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under BHP limit with current reservoir condition
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under THP limit with current reservoir condition
void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
// update WellState based on IPR and associated VFP table
void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
// for a well, when all drawdown are in the wrong direction, then this well will not
// be able to produce/inject .
bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
// whether the well can produce / inject based on the current well state (bhp)
bool canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger);
// turn on crossflow to avoid singular well equations
// when the well is banned from cross-flow and the BHP is not properly initialized,
// we turn on crossflow to avoid singular well equations. It can result in wrong-signed
// well rates, it can cause problem for THP calculation
// TODO: looking for better alternative to avoid wrong-signed well rates
bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
// calculate the BHP from THP target based on IPR
// TODO: we need to check the operablility here first, if not operable, then maybe there is
// no point to do this
double calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const;
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
const double dx);
// calculate a relaxation factor to avoid overshoot of the fractions for producers
// which might result in negative rates
static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override;
// calculate the skin pressure based on water velocity, throughput and polymer concentration.
// throughput is used to describe the formation damage during water/polymer injection.
// calculated skin pressure will be applied to the drawdown during perforation rate calculation
// to handle the effect from formation damage.
EvalWell pskin(const double throuhgput,
const EvalWell& water_velocity,
const EvalWell& poly_inj_conc,
Opm::DeferredLogger& deferred_logger) const;
// calculate the skin pressure based on water velocity, throughput during water injection.
EvalWell pskinwater(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const;
// calculate the injecting polymer molecular weight based on the througput and water velocity
EvalWell wpolymermw(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const;
// handle the extra equations for polymer injectivity study
void handleInjectivityRateAndEquations(const IntensiveQuantities& int_quants,
const WellState& well_state,
const int perf,
std::vector<EvalWell>& cq_s,
Opm::DeferredLogger& deferred_logger);
virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
void checkConvergenceControlEq(ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
};
}
#include "StandardWellV_impl.hpp"
#endif // OPM_STANDARDWELLV_HEADER_INCLUDED

3182
opm/simulators/wells/StandardWellV_impl.hpp
View File

File diff suppressed because it is too large Load Diff

426
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -34,8 +34,6 @@ namespace Opm
: Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components)
, perf_densities_(number_of_perforations_)
, perf_pressure_diffs_(number_of_perforations_)
, primary_variables_(numWellEq, 0.0)
, primary_variables_evaluation_(numWellEq) // the number of the primary variables
, F0_(numWellConservationEq)
, ipr_a_(number_of_phases_)
, ipr_b_(number_of_phases_)
@ -67,6 +65,18 @@ namespace Opm
perf_depth_[perf] = depth_arg[cell_idx];
}
// counting/updating primary variable numbers
if (this->has_polymermw && well_type_ == INJECTOR) {
// adding a primary variable for water perforation rate per connection
numWellEq_ += number_of_perforations_;
// adding a primary variable for skin pressure per connection
numWellEq_ += number_of_perforations_;
}
// with the updated numWellEq_, we can initialize the primary variables and matrices now
primary_variables_.resize(numWellEq_, 0.0);
primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + numEq, 0.0});
// setup sparsity pattern for the matrices
//[A C^T [x = [ res
// B D] x_well] res_well]
@ -79,6 +89,8 @@ namespace Opm
// Add nonzeros for diagonal
row.insert(row.index());
}
// the block size is run-time determined now
invDuneD_[0][0].resize(numWellEq_, numWellEq_);
for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
for (int perf = 0 ; perf < number_of_perforations_; ++perf) {
@ -87,6 +99,12 @@ namespace Opm
}
}
for (int perf = 0 ; perf < number_of_perforations_; ++perf) {
const int cell_idx = well_cells_[perf];
// the block size is run-time determined now
duneB_[0][cell_idx].resize(numWellEq_, numEq);
}
// make the C^T matrix
for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
@ -95,11 +113,25 @@ namespace Opm
}
}
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int cell_idx = well_cells_[perf];
duneC_[0][cell_idx].resize(numWellEq_, numEq);
}
resWell_.resize(1);
// the block size of resWell_ is also run-time determined now
resWell_[0].resize(numWellEq_);
// resize temporary class variables
Bx_.resize( duneB_.N() );
for (unsigned i = 0; i < duneB_.N(); ++i) {
Bx_[i].resize(numWellEq_);
}
invDrw_.resize( invDuneD_.N() );
for (unsigned i = 0; i < invDuneD_.N(); ++i) {
invDrw_[i].resize(numWellEq_);
}
}
@ -110,12 +142,9 @@ namespace Opm
void StandardWell<TypeTag>::
initPrimaryVariablesEvaluation() const
{
for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
assert( (size_t)eqIdx < primary_variables_.size() );
primary_variables_evaluation_[eqIdx] = 0.0;
primary_variables_evaluation_[eqIdx].setValue(primary_variables_[eqIdx]);
primary_variables_evaluation_[eqIdx].setDerivative(numEq + eqIdx, 1.0);
for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) {
primary_variables_evaluation_[eqIdx] =
EvalWell::createVariable(numWellEq_ + numEq, primary_variables_[eqIdx], numEq + eqIdx);
}
}
@ -221,7 +250,7 @@ namespace Opm
}
// Oil fraction
EvalWell well_fraction = 1.0;
EvalWell well_fraction(numWellEq_ + numEq, 1.0);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
well_fraction -= primary_variables_evaluation_[WFrac];
}
@ -245,7 +274,7 @@ namespace Opm
wellSurfaceVolumeFraction(const int compIdx) const
{
EvalWell sum_volume_fraction_scaled = 0.;
EvalWell sum_volume_fraction_scaled(numWellEq_ + numEq, 0.);
for (int idx = 0; idx < num_components_; ++idx) {
sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
}
@ -264,8 +293,7 @@ namespace Opm
StandardWell<TypeTag>::
extendEval(const Eval& in) const
{
EvalWell out = 0.0;
out.setValue(in.value());
EvalWell out(numWellEq_ + numEq, in.value());
for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
out.setDerivative(eqIdx, in.derivative(eqIdx));
}
@ -290,11 +318,12 @@ namespace Opm
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const
{
const auto& fs = intQuants.fluidState();
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(num_components_, 0.0);
std::vector<EvalWell> b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0});
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@ -309,7 +338,13 @@ namespace Opm
// Pressure drawdown (also used to determine direction of flow)
const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
const EvalWell drawdown = pressure - well_pressure;
EvalWell drawdown = pressure - well_pressure;
if (this->has_polymermw && well_type_ == INJECTOR) {
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
drawdown += skin_pressure;
}
// producing perforations
if ( drawdown.value() > 0 ) {
@ -358,13 +393,13 @@ namespace Opm
const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
// surface volume fraction of fluids within wellbore
std::vector<EvalWell> cmix_s(num_components_, 0.);
std::vector<EvalWell> cmix_s(num_components_, EvalWell{numWellEq_ + numEq});
for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) {
cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
}
// compute volume ratio between connection at standard conditions
EvalWell volumeRatio = 0.0;
EvalWell volumeRatio(numWellEq_ + numEq, 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
@ -378,7 +413,7 @@ namespace Opm
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RS/RV factors if both oil and gas active
const EvalWell d = 1.0 - rv * rs;
const EvalWell d = EvalWell(numWellEq_ + numEq, 1.0) - rv * rs;
if (d.value() == 0.0) {
OPM_DEFLOG_THROW(Opm::NumericalIssue, "Zero d value obtained for well " << name() << " during flux calcuation"
@ -449,7 +484,11 @@ namespace Opm
WellState& well_state,
Opm::DeferredLogger& deferred_logger)
{
<<<<<<< 6cd9751a394b4d8f839c1b3b7201e4ddc2bcf664:opm/simulators/wells/StandardWell_impl.hpp
SummaryState summaryState;
=======
// TODO: only_wells should be put back to save some computation
>>>>>>> removing StandardWellV:opm/autodiff/StandardWell_impl.hpp
checkWellOperability(ebosSimulator, well_state, deferred_logger);
@ -481,10 +520,10 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(num_components_, 0.0);
std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.});
getMobility(ebosSimulator, perf, mob, deferred_logger);
std::vector<EvalWell> cq_s(num_components_, 0.0);
std::vector<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
@ -492,6 +531,11 @@ namespace Opm
computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
// better way to do here is that use the cq_s and then replace the cq_s_water here?
if (has_polymer && this->has_polymermw && well_type_ == INJECTOR) {
handleInjectivityRateAndEquations(intQuants, well_state, perf, cq_s, deferred_logger);
}
// updating the solution gas rate and solution oil rate
if (well_type_ == PRODUCER) {
well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
@ -512,7 +556,7 @@ namespace Opm
resWell_[0][componentIdx] += cq_s_effective.value();
// assemble the jacobians
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
// also need to consider the efficiency factor when manipulating the jacobians.
duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
invDuneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+numEq);
@ -539,7 +583,7 @@ namespace Opm
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
// convert to reservoar conditions
EvalWell cq_r_thermal = 0.0;
EvalWell cq_r_thermal(numWellEq_ + numEq, 0.);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if(FluidSystem::waterPhaseIdx == phaseIdx)
@ -601,8 +645,21 @@ namespace Opm
connectionRates_[perf][contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
if (this->has_polymermw) {
if (well_type_ == PRODUCER) {
// the source term related to transport of molecular weight
EvalWell cq_s_polymw = cq_s_poly;
if (well_type_ == INJECTOR) {
const int wat_vel_index = Bhp + 1 + perf;
const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index];
if (water_velocity > 0.) { // injecting
const double throughput = well_state.perfThroughput()[first_perf_ + perf];
const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
cq_s_polymw *= molecular_weight;
} else {
// we do not consider the molecular weight from the polymer
// going-back to the wellbore through injector
cq_s_polymw *= 0.;
}
} else if (well_type_ == PRODUCER) {
if (cq_s_polymw < 0.) {
cq_s_polymw *= extendEval(intQuants.polymerMoleWeight() );
} else {
@ -610,8 +667,8 @@ namespace Opm
// re-injecting back through producer
cq_s_polymw *= 0.;
}
connectionRates_[perf][this->contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
}
connectionRates_[perf][this->contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
}
}
@ -645,7 +702,7 @@ namespace Opm
// since all the rates are under surface condition
EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
resWell_loc -= getQs(componentIdx) * well_efficiency_factor_;
for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
}
resWell_[0][componentIdx] += resWell_loc.value();
@ -672,11 +729,11 @@ namespace Opm
StandardWell<TypeTag>::
assembleControlEq(Opm::DeferredLogger& deferred_logger)
{
EvalWell control_eq(0.0);
EvalWell control_eq(numWellEq_ + numEq, 0.);
switch (well_controls_get_current_type(well_controls_)) {
case THP:
{
std::vector<EvalWell> rates(3, 0.);
std::vector<EvalWell> rates(3, {numWellEq_ + numEq, 0.});
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = getQs(flowPhaseToEbosCompIdx(Water));
}
@ -703,7 +760,7 @@ namespace Opm
control_eq = getWQTotal() - target_rate;
} else if (well_type_ == PRODUCER) {
if (target_rate != 0.) {
EvalWell rate_for_control(0.);
EvalWell rate_for_control(numWellEq_ + numEq, 0.);
const EvalWell& g_total = getWQTotal();
// a variable to check if we are producing any targeting fluids
double sum_fraction = 0.;
@ -752,7 +809,7 @@ namespace Opm
}
} else {
const EvalWell& g_total = getWQTotal();
EvalWell rate_for_control(0.0); // reservoir rate
EvalWell rate_for_control(numWellEq_ + numEq, 0.); // reservoir rate
for (int phase = 0; phase < number_of_phases_; ++phase) {
rate_for_control += g_total * wellVolumeFraction( flowPhaseToEbosCompIdx(phase) );
}
@ -767,7 +824,7 @@ namespace Opm
// using control_eq to update the matrix and residuals
// TODO: we should use a different index system for the well equations
resWell_[0][Bhp] = control_eq.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) {
invDuneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + numEq);
}
}
@ -913,6 +970,21 @@ namespace Opm
// 1e5 to make sure bhp will not be below 1bar
primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
}
// for the water velocity and skin pressure
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const double relaxation_factor = 0.9;
const double dx_wat_vel = dwells[0][wat_vel_index];
primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
const double dx_pskin = dwells[0][pskin_index];
primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
}
}
}
@ -1067,6 +1139,14 @@ namespace Opm
}
updateThp(well_state, deferred_logger);
// other primary variables related to polymer injectivity study
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
}
}
}
@ -1121,7 +1201,6 @@ namespace Opm
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
{
// number of phases
const int np = number_of_phases_;
const int well_index = index_of_well_;
@ -1241,7 +1320,7 @@ namespace Opm
std::fill(ipr_b_.begin(), ipr_b_.end(), 0.);
for (int perf = 0; perf < number_of_perforations_; ++perf) {
std::vector<EvalWell> mob(num_components_, 0.0);
std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.0});
// TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobility(ebos_simulator, perf, mob, deferred_logger);
@ -1324,8 +1403,7 @@ namespace Opm
StandardWell<TypeTag>::
checkWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
)
Opm::DeferredLogger& deferred_logger)
{
// focusing on PRODUCER for now
if (well_type_ == INJECTOR) {
@ -1358,8 +1436,7 @@ namespace Opm
StandardWell<TypeTag>::
updateWellOperability(const Simulator& ebos_simulator,
const WellState& /* well_state */,
Opm::DeferredLogger& deferred_logger
)
Opm::DeferredLogger& deferred_logger)
{
this->operability_status_.reset();
@ -1405,7 +1482,7 @@ namespace Opm
// option 2: stick with the above IPR curve
// we use IPR here
std::vector<double> well_rates_bhp_limit;
computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp_limit), well_rates_bhp_limit, deferred_logger);
const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger);
const double thp_limit = this->getTHPConstraint(deferred_logger);
@ -1506,8 +1583,7 @@ namespace Opm
StandardWell<TypeTag>::
canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
)
Opm::DeferredLogger& deferred_logger)
{
const double bhp = well_state.bhp()[index_of_well_];
std::vector<double> well_rates;
@ -1595,7 +1671,7 @@ namespace Opm
initPrimaryVariablesEvaluation();
std::vector<double> rates;
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp), rates, deferred_logger);
// TODO: double checke the obtained rates
// this is another places we might obtain negative rates
@ -1920,8 +1996,8 @@ namespace Opm
const double tol_wells = param_.tolerance_wells_;
const double maxResidualAllowed = param_.max_residual_allowed_;
std::vector<double> res(numWellEq);
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
std::vector<double> res(numWellEq_);
for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
// magnitude of the residual matters
res[eq_idx] = std::abs(resWell_[0][eq_idx]);
}
@ -1956,7 +2032,70 @@ namespace Opm
}
}
<<<<<<< 6cd9751a394b4d8f839c1b3b7201e4ddc2bcf664:opm/simulators/wells/StandardWell_impl.hpp
checkConvergenceControlEq(report, deferred_logger);
=======
// processing the residual of the well control equation
const double well_control_residual = res[numWellEq_ - 1];
// TODO: we should have better way to specify the control equation tolerance
double control_tolerance = 0.;
switch(well_controls_get_current_type(well_controls_)) {
case THP:
type = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case BHP: // pressure type of control
type = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case RESERVOIR_RATE:
case SURFACE_RATE:
type = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; // smaller tolerance for rate control
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " + name(), deferred_logger);
}
const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({type, CR::Severity::NotANumber, dummy_component, name()});
} else if (well_control_residual > maxResidualAllowed * 10.) {
report.setWellFailed({type, CR::Severity::TooLarge, dummy_component, name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
}
>>>>>>> removing StandardWellV:opm/autodiff/StandardWell_impl.hpp
if (this->has_polymermw && well_type_ == INJECTOR) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, name()});
}
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double pskin_residual = res[Bhp + 1 + perf + number_of_perforations_];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, name()});
}
}
}
return report;
}
@ -2029,6 +2168,7 @@ namespace Opm
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
BVectorWell dx_well(1);
dx_well[0].resize(numWellEq_);
invDuneD_.mv(resWell_, dx_well);
updateWellState(dx_well, well_state, deferred_logger);
@ -2146,6 +2286,8 @@ namespace Opm
if (!this->isOperable()) return;
BVectorWell xw(1);
xw[0].resize(numWellEq_);
recoverSolutionWell(x, xw);
updateWellState(xw, well_state, deferred_logger);
}
@ -2171,12 +2313,12 @@ namespace Opm
const int cell_idx = well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
// flux for each perforation
std::vector<EvalWell> mob(num_components_, 0.0);
std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.});
getMobility(ebosSimulator, perf, mob, deferred_logger);
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = well_index_[perf] * trans_mult;
std::vector<EvalWell> cq_s(num_components_, 0.0);
std::vector<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, EvalWell(bhp), Tw, perf, allow_cf,
@ -2308,7 +2450,11 @@ namespace Opm
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
<<<<<<< 6cd9751a394b4d8f839c1b3b7201e4ddc2bcf664:opm/simulators/wells/StandardWell_impl.hpp
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, potentials, deferred_logger);
=======
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), potentials, deferred_logger);
>>>>>>> removing StandardWellV:opm/autodiff/StandardWell_impl.hpp
// checking whether the potentials have valid values
for (const double value : potentials) {
@ -2355,7 +2501,6 @@ namespace Opm
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) // const
{
updatePrimaryVariables(well_state, deferred_logger);
computeWellConnectionPressures(ebosSimulator, well_state);
@ -2372,7 +2517,12 @@ namespace Opm
// does the well have a THP related constraint?
if ( !wellHasTHPConstraints() ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
<<<<<<< 6cd9751a394b4d8f839c1b3b7201e4ddc2bcf664:opm/simulators/wells/StandardWell_impl.hpp
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
=======
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger);
>>>>>>> removing StandardWellV:opm/autodiff/StandardWell_impl.hpp
} else {
// the well has a THP related constraint
// checking whether a well is newly added, it only happens at the beginning of the report step
@ -2384,7 +2534,11 @@ namespace Opm
}
} else {
// We need to generate a reasonable rates to start the iteration process
<<<<<<< 6cd9751a394b4d8f839c1b3b7201e4ddc2bcf664:opm/simulators/wells/StandardWell_impl.hpp
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
=======
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger);
>>>>>>> removing StandardWellV:opm/autodiff/StandardWell_impl.hpp
for (double& value : well_potentials) {
// make the value a little safer in case the BHP limits are default ones
// TODO: a better way should be a better rescaling based on the investigation of the VFP table.
@ -2406,7 +2560,6 @@ namespace Opm
StandardWell<TypeTag>::
updatePrimaryVariables(const WellState& well_state, Opm::DeferredLogger& deferred_logger) const
{
if (!this->isOperable()) return;
const int well_index = index_of_well_;
@ -2489,6 +2642,14 @@ namespace Opm
// BHP
primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
// other primary variables related to polymer injectio
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
}
}
}
@ -2523,26 +2684,23 @@ namespace Opm
// TODO: it is possible it should be a Evaluation
const double rho = perf_densities_[0];
ValueType bhp = 0.;
if (well_type_ == INJECTOR) {
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
return vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
}
else if (well_type_ == PRODUCER) {
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
return vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
return bhp;
}
@ -2604,6 +2762,11 @@ namespace Opm
std::vector<EvalWell>& mob,
Opm::DeferredLogger& deferred_logger) const
{
// for the cases related to polymer molecular weight, we assume fully mixing
// as a result, the polymer and water share the same viscosity
if (this->has_polymermw) {
return;
}
const int cell_idx = well_cells_[perf];
const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const EvalWell polymer_concentration = extendEval(int_quant.polymerConcentration());
@ -2628,7 +2791,7 @@ namespace Opm
const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
const EvalWell& bhp = getBhp();
std::vector<EvalWell> cq_s(num_components_,0.0);
std::vector<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
@ -2684,9 +2847,9 @@ namespace Opm
while ( col != row.end() && col.index() < col_index ) ++col;
assert(col != row.end() && col.index() == col_index);
Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp;
Dune::DynamicMatrix<Scalar> tmp;
Detail::multMatrix(invDuneD_[0][0], (*colB), tmp);
typename Mat::block_type tmp1;
Dune::FMatrixHelp::multMatrix(invDuneD_[0][0], (*colB), tmp);
Detail::multMatrixTransposed((*colC), tmp, tmp1);
(*col) -= tmp1;
}
@ -2865,11 +3028,166 @@ namespace Opm
template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>::
pskinwater(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting skin pressure is requested for well " << name(), deferred_logger);
}
const int water_table_id = well_ecl_->getPolymerProperties(current_step_).m_skprwattable;
if (water_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
}
const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_water(numWellEq_ + numEq, 0.0);
pskin_water = water_table_func.eval(throughput_eval, water_velocity);
return pskin_water;
}
template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>::
pskin(const double throughput,
const EvalWell& water_velocity,
const EvalWell& poly_inj_conc,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting skin pressure is requested for well " << name(), deferred_logger);
}
const double sign = water_velocity >= 0. ? 1.0 : -1.0;
const EvalWell water_velocity_abs = Opm::abs(water_velocity);
if (poly_inj_conc == 0.) {
return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
}
const int polymer_table_id = well_ecl_->getPolymerProperties(current_step_).m_skprpolytable;
if (polymer_table_id <= 0) {
OPM_DEFLOG_THROW(std::runtime_error, "Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
}
const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
const double reference_concentration = skprpolytable.refConcentration;
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
// the skin pressure when injecting water, which also means the polymer concentration is zero
EvalWell pskin_poly(numWellEq_ + numEq, 0.0);
pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
if (poly_inj_conc == reference_concentration) {
return sign * pskin_poly;
}
// poly_inj_conc != reference concentration of the table, then some interpolation will be required
const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
return sign * pskin;
}
template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>::
wpolymermw(const double throughput,
const EvalWell& water_velocity,
Opm::DeferredLogger& deferred_logger) const
{
if (!this->has_polymermw) {
OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
"while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
}
const int table_id = well_ecl_->getPolymerProperties(current_step_).m_plymwinjtable;
const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
const EvalWell throughput_eval(numWellEq_ + numEq, throughput);
EvalWell molecular_weight(numWellEq_ + numEq, 0.);
if (wpolymer() == 0.) { // not injecting polymer
return molecular_weight;
}
molecular_weight = table_func.eval(throughput_eval, Opm::abs(water_velocity));
return molecular_weight;
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
updateWaterThroughput(const double dt OPM_UNUSED, WellState& well_state OPM_UNUSED) const
updateWaterThroughput(const double dt, WellState &well_state) const
{
if (this->has_polymermw && well_type_ == INJECTOR) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
if (perf_water_vel > 0.) {
well_state.perfThroughput()[first_perf_ + perf] += perf_water_vel * dt;
}
}
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
handleInjectivityRateAndEquations(const IntensiveQuantities& int_quants,
const WellState& well_state,
const int perf,
std::vector<EvalWell>& cq_s,
Opm::DeferredLogger& deferred_logger)
{
const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const EvalWell& water_flux_s = cq_s[water_comp_idx];
const auto& fs = int_quants.fluidState();
const EvalWell b_w = extendEval(fs.invB(FluidSystem::waterPhaseIdx));
const EvalWell water_flux_r = water_flux_s / b_w;
const double area = M_PI * bore_diameters_[perf] * perf_length_[perf];
const EvalWell water_velocity = water_flux_r / area;
const int wat_vel_index = Bhp + 1 + perf;
// equation for the water velocity
const EvalWell eq_wat_vel = primary_variables_evaluation_[wat_vel_index] - water_velocity;
resWell_[0][wat_vel_index] = eq_wat_vel.value();
const double throughput = well_state.perfThroughput()[first_perf_ + perf];
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
EvalWell poly_conc(numWellEq_ + numEq, 0.0);
poly_conc.setValue(wpolymer());
// equation for the skin pressure
const EvalWell eq_pskin = primary_variables_evaluation_[pskin_index]
- pskin(throughput, primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
resWell_[0][pskin_index] = eq_pskin.value();
for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
invDuneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+numEq);
invDuneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+numEq);
}
// water rate is update to use the form from water velocity, since water velocity is
// a primary variable now
cq_s[water_comp_idx] = area * primary_variables_evaluation_[wat_vel_index] * b_w;
// the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
const int cell_idx = well_cells_[perf];
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
}
}

4
tests/test_wellmodel.cpp
View File

@ -208,7 +208,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
BOOST_CHECK_EQUAL(well->name(), "PROD1");
BOOST_CHECK(well->wellType() == PRODUCER);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numWellEq == 4);
BOOST_CHECK(well->numStaticWellEq== 4);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);
@ -227,7 +227,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->wellType() == INJECTOR);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numWellEq == 4);
BOOST_CHECK(well->numStaticWellEq== 4);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);