OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1250 from GitPaean/refactoring_well_model_June_2017_rebase

Browse Source 
well refactoring for flow_ebos
This commit is contained in:
Atgeirr Flø Rasmussen
2017-08-28 14:36:18 +02:00
committed by GitHub
parent 060631c5dc 9accb56c86
commit 160314424e
15 changed files with 4166 additions and 3078 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists_files.cmake
View File

@@ -91,6 +91,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_solventprops_ad.cpp
tests/test_multisegmentwells.cpp
tests/test_multiphaseupwind.cpp
tests/test_wellmodel.cpp
# tests/test_thresholdpressure.cpp
tests/test_wellswitchlogger.cpp
tests/test_timer.cpp
@@ -102,6 +103,7 @@ list (APPEND TEST_DATA_FILES
tests/VFPPROD2
tests/msw.data
tests/TESTTIMER.DATA
tests/TESTWELLMODEL.DATA
)
@@ -239,6 +241,10 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/WellHelpers.hpp
opm/autodiff/StandardWells.hpp
opm/autodiff/StandardWells_impl.hpp
opm/autodiff/WellInterface.hpp
opm/autodiff/WellInterface_impl.hpp
opm/autodiff/StandardWell.hpp
opm/autodiff/StandardWell_impl.hpp
opm/autodiff/StandardWellsDense.hpp
opm/autodiff/StandardWellsSolvent.hpp
opm/autodiff/StandardWellsSolvent_impl.hpp

76
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -154,7 +154,7 @@ namespace Opm {
/// \param[in] terminal_output request output to cout/cerr
BlackoilModelEbos(Simulator& ebosSimulator,
const ModelParameters& param,
const StandardWellsDense<TypeTag>& well_model,
StandardWellsDense<TypeTag>& well_model,
RateConverterType& rate_converter,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output
@@ -284,10 +284,9 @@ namespace Opm {
const int nc = AutoDiffGrid::numCells(grid_);
const int nw = numWells();
BVector x(nc);
BVector xw(nw);
try {
solveJacobianSystem(x, xw);
solveJacobianSystem(x);
report.linear_solve_time += perfTimer.stop();
report.total_linear_iterations += linearIterationsLastSolve();
}
@@ -322,7 +321,11 @@ namespace Opm {
// Apply the update, with considering model-dependent limitations and
// chopping of the update.
updateState(x,iteration);
wellModel().updateWellState(xw, well_state);
if( nw > 0 )
{
wellModel().recoverWellSolutionAndUpdateWellState(x, well_state);
}
report.update_time += perfTimer.stop();
}
@@ -488,7 +491,7 @@ namespace Opm {
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
void solveJacobianSystem(BVector& x, BVector& xw) const
void solveJacobianSystem(BVector& x) const
{
const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
@@ -510,13 +513,6 @@ namespace Opm {
Operator opA(ebosJac, well_model_);
istlSolver().solve( opA, x, ebosResid );
}
if( xw.size() > 0 )
{
// recover wells.
xw = 0.0;
wellModel().recoverVariable(x, xw);
}
}
//=====================================================================
@@ -766,8 +762,7 @@ namespace Opm {
const double pvSumLocal,
std::vector< Scalar >& R_sum,
std::vector< Scalar >& maxCoeff,
std::vector< Scalar >& B_avg,
std::vector< Scalar >& maxNormWell )
std::vector< Scalar >& B_avg)
{
// Compute total pore volume (use only owned entries)
double pvSum = pvSumLocal;
@@ -779,13 +774,12 @@ namespace Opm {
std::vector< Scalar > maxBuffer;
const int numComp = B_avg.size();
sumBuffer.reserve( 2*numComp + 1 ); // +1 for pvSum
maxBuffer.reserve( 2*numComp );
maxBuffer.reserve( numComp );
for( int compIdx = 0; compIdx < numComp; ++compIdx )
{
sumBuffer.push_back( B_avg[ compIdx ] );
sumBuffer.push_back( R_sum[ compIdx ] );
maxBuffer.push_back( maxCoeff[ compIdx ] );
maxBuffer.push_back( maxNormWell[ compIdx ] );
}
// Compute total pore volume
@@ -801,11 +795,14 @@ namespace Opm {
for( int compIdx = 0, buffIdx = 0; compIdx < numComp; ++compIdx, ++buffIdx )
{
B_avg[ compIdx ] = sumBuffer[ buffIdx ];
maxCoeff[ compIdx ] = maxBuffer[ buffIdx ];
++buffIdx;
R_sum[ compIdx ] = sumBuffer[ buffIdx ];
maxNormWell[ compIdx ] = maxBuffer[ buffIdx ];
}
for( int compIdx = 0; compIdx < numComp; ++compIdx )
{
maxCoeff[ compIdx ] = maxBuffer[ compIdx ];
}
// restore global pore volume
@@ -828,7 +825,6 @@ namespace Opm {
const double dt = timer.currentStepLength();
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
const double tol_wells = param_.tolerance_wells_;
const int np = numPhases();
const int numComp = numComponents();
@@ -836,7 +832,6 @@ namespace Opm {
Vector R_sum(numComp, 0.0 );
Vector B_avg(numComp, 0.0 );
Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
Vector maxNormWell(numComp, 0.0 );
const auto& ebosModel = ebosSimulator_.model();
const auto& ebosProblem = ebosSimulator_.problem();
@@ -896,24 +891,16 @@ namespace Opm {
B_avg[ i ] /= Scalar( global_nc_ );
}
// compute maximum of local well residuals
const Vector& wellResidual = wellModel().residual();
const int nw = wellResidual.size() / numComp;
assert(nw * numComp == int(wellResidual.size()));
for( int compIdx = 0; compIdx < numComp; ++compIdx )
{
for ( int w = 0; w < nw; ++w ) {
maxNormWell[compIdx] = std::max(maxNormWell[compIdx], std::abs(wellResidual[nw*compIdx + w]));
}
}
// TODO: we remove the maxNormWell for now because the convergence of wells are on a individual well basis.
// Anyway, we need to provide some infromation to help debug the well iteration process.
// compute global sum and max of quantities
const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
R_sum, maxCoeff, B_avg, maxNormWell );
R_sum, maxCoeff, B_avg);
Vector CNV(numComp);
Vector mass_balance_residual(numComp);
Vector well_flux_residual(numComp);
bool converged_MB = true;
bool converged_CNV = true;
@@ -927,8 +914,7 @@ namespace Opm {
converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
// Well flux convergence is only for fluid phases, not other materials
// in our current implementation.
well_flux_residual[compIdx] = B_avg[compIdx] * maxNormWell[compIdx];
converged_Well = converged_Well && (well_flux_residual[compIdx] < tol_wells);
converged_Well = wellModel().getWellConvergence(ebosSimulator_, B_avg);
residual_norms.push_back(CNV[compIdx]);
}
@@ -965,9 +951,6 @@ namespace Opm {
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " CNV(" + key[ compIdx ] + ") ";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " W-FLUX(" + key[ compIdx ] + ")";
}
OpmLog::note(msg);
}
std::ostringstream ss;
@@ -980,9 +963,6 @@ namespace Opm {
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << CNV[compIdx];
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
ss << std::setw(11) << well_flux_residual[compIdx];
}
ss.precision(oprec);
ss.flags(oflags);
OpmLog::note(ss.str());
@@ -992,13 +972,11 @@ namespace Opm {
const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
if (std::isnan(mass_balance_residual[phaseIdx])
|| std::isnan(CNV[phaseIdx])
|| (phaseIdx < numPhases() && std::isnan(well_flux_residual[phaseIdx]))) {
|| std::isnan(CNV[phaseIdx])) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName);
}
if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
|| CNV[phaseIdx] > maxResidualAllowed()
|| (phaseIdx < numPhases() && well_flux_residual[phaseIdx] > maxResidualAllowed())) {
|| CNV[phaseIdx] > maxResidualAllowed()) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName);
}
}
@@ -1523,7 +1501,7 @@ namespace Opm {
SimulatorReport failureReport_;
// Well Model
StandardWellsDense<TypeTag> well_model_;
StandardWellsDense<TypeTag>& well_model_;
/// \brief Whether we print something to std::cout
bool terminal_output_;
@@ -1545,13 +1523,7 @@ namespace Opm {
const StandardWellsDense<TypeTag>&
wellModel() const { return well_model_; }
/// return the Well struct in the StandardWells
const Wells& wells() const { return well_model_.wells(); }
/// return true if wells are available in the reservoir
bool wellsActive() const { return well_model_.wellsActive(); }
int numWells() const { return wellsActive() ? wells().number_of_wells : 0; }
int numWells() const { return well_model_.numWells(); }
/// return true if wells are available on this process
bool localWellsActive() const { return well_model_.localWellsActive(); }

3
opm/autodiff/Compat.cpp
View File

@@ -294,9 +294,6 @@ void wellsToState( const data::Wells& wells,
{
// Set base class variables.
wellsToState(wells, phases, static_cast<WellStateFullyImplicitBlackoil&>(state));
// Set wellSolution() variable.
state.setWellSolutions(phases);
}

19
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -141,7 +141,6 @@ public:
ExtraData extra;
failureReport_ = SimulatorReport();
extractLegacyPoreVolume_();
extractLegacyDepth_();
// communicate the initial solution to ebos
@@ -480,7 +479,6 @@ protected:
activePhases,
gravity,
legacyDepth_,
legacyPoreVolume_,
globalNumCells,
grid());
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
@@ -867,22 +865,6 @@ protected:
}
}
void extractLegacyPoreVolume_()
{
const auto& grid = ebosSimulator_.gridManager().grid();
const unsigned numCells = grid.size(/*codim=*/0);
const auto& ebosProblem = ebosSimulator_.problem();
const auto& ebosModel = ebosSimulator_.model();
legacyPoreVolume_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
// todo (?): respect rock compressibility
legacyPoreVolume_[cellIdx] =
ebosModel.dofTotalVolume(cellIdx)
*ebosProblem.porosity(cellIdx);
}
}
void extractLegacyDepth_()
{
const auto& grid = ebosSimulator_.gridManager().grid();
@@ -1009,7 +991,6 @@ protected:
Simulator& ebosSimulator_;
std::vector<int> legacyCellPvtRegionIdx_;
std::vector<double> legacyPoreVolume_;
std::vector<double> legacyDepth_;
typedef typename Solver::SolverParameters SolverParameters;

305
opm/autodiff/StandardWell.hpp Normal file
View File

@@ -0,0 +1,305 @@
/*
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_STANDARDWELL_HEADER_INCLUDED
#define OPM_STANDARDWELL_HEADER_INCLUDED
#include <opm/autodiff/WellInterface.hpp>
namespace Opm
{
template<typename TypeTag>
class StandardWell: 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::BlackoilIndices;
using typename Base::PolymerModule;
// the positions of the primary variables for StandardWell
// there are three primary variables, the second and the third ones are F_w and F_g
// the first one can be total rate (G_t) or bhp, based on the control
enum WellVariablePositions {
XvarWell = 0,
WFrac = 1,
GFrac = 2,
SFrac = 3
};
using typename Base::Scalar;
using typename Base::ConvergenceReport;
using Base::numEq;
using Base::has_solvent;
using Base::has_polymer;
using Base::name;
// TODO: with flow_ebos，for a 2P deck, // TODO: for the 2p deck, numEq will be 3, a dummy phase is already added from the reservoir side.
// it will cause problem here without processing the dummy phase.
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // number of wellEq is only numEq - 1 for polymer
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;
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
// TODO: should these go to WellInterface?
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
StandardWell(const Well* well, const int time_step, const Wells* wells);
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& 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);
/// 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, const ModelParameters& param,
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;
virtual void updatePrimaryVariables(const WellState& well_state) const;
protected:
// protected functions from the Base class
using Base::getAllowCrossFlow;
using Base::phaseUsage;
using Base::active;
using Base::flowToEbosPvIdx;
using Base::flowPhaseToEbosPhaseIdx;
using Base::flowPhaseToEbosCompIdx;
using Base::numComponents;
using Base::wsolvent;
using Base::wpolymer;
using Base::wellHasTHPConstraints;
using Base::mostStrictBhpFromBhpLimits;
// protected member variables from the Base class
using Base::vfp_properties_;
using Base::gravity_;
using Base::well_efficiency_factor_;
using Base::phase_usage_;
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::perf_rep_radius_;
using Base::perf_length_;
using Base::bore_diameters_;
// 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_;
// TODO: this function should be moved to the base class.
// while it faces chanllenges for MSWell later, since the calculation of bhp
// based on THP is never implemented for MSWell yet.
EvalWell getBhp() const;
// TODO: it is also possible to be moved to the base class.
EvalWell getQs(const int comp_idx) const;
EvalWell wellVolumeFractionScaled(const int phase) const;
EvalWell wellVolumeFraction(const int phase) const;
EvalWell wellSurfaceVolumeFraction(const int phase) const;
EvalWell extendEval(const Eval& in) const;
bool crossFlowAllowed(const Simulator& ebosSimulator) 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,
const BlackoilModelParameters& param,
WellState& well_state) const;
// calculate the properties for the well connections
// to calulate the pressure difference between well connections.
void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw,
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& xw,
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);
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state);
// TODO: to check whether all the paramters are required
void computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob_perfcells_dense,
const double Tw, const EvalWell& bhp, const double& cdp,
const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
// TODO: maybe we should provide a light version of computePerfRate, which does not include the
// calculation of the derivatives
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
std::vector<double>& well_flux) const;
std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential) const;
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index) const;
double calculateThpFromBhp(const std::vector<double>& rates, const int control_index, const double bhp) const;
// get the mobility for specific perforation
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob) const;
};
}
#include "StandardWell_impl.hpp"
#endif // OPM_STANDARDWELL_HEADER_INCLUDED

2060
opm/autodiff/StandardWell_impl.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

438
opm/autodiff/StandardWellsDense.hpp
View File

@@ -39,8 +39,6 @@
#include <opm/core/wells/WellCollection.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
@@ -49,27 +47,19 @@
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/WellInterface.hpp>
#include <opm/autodiff/StandardWell.hpp>
#include<dune/common/fmatrix.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/matrixmatrix.hh>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
#include <math.h>
namespace Opm {
enum WellVariablePositions {
XvarWell = 0,
WFrac = 1,
GFrac = 2,
SFrac = 3
};
/// Class for handling the standard well model.
template<typename TypeTag>
class StandardWellsDense {
@@ -82,26 +72,17 @@ enum WellVariablePositions {
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
static const int numEq = BlackoilIndices::numEq;
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // //numEq; //number of wellEq is only numEq for polymer
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
// TODO: where we should put these types, WellInterface or Well Model?
// or there is some other strategy, like TypeTag
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<Scalar, /*size=*/numEq + numWellEq> EvalWell;
typedef DenseAd::Evaluation<Scalar, /*size=*/numEq> Eval;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
// For the conversion between the surface volume rate and resrevoir voidage rate
@@ -121,13 +102,143 @@ enum WellVariablePositions {
const std::vector<bool>& active_arg,
const double gravity_arg,
const std::vector<double>& depth_arg,
const std::vector<double>& pv_arg,
long int global_nc,
const Grid& grid);
void setVFPProperties(const VFPProperties* vfp_properties_arg);
/// The number of components in the model.
SimulatorReport assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state);
// substract Binv(D)rw from r;
void apply( BVector& r) const;
// subtract B*inv(D)*C * x from A*x
void apply(const BVector& x, BVector& Ax) const;
// apply well model with scaling of alpha
void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
// using the solution x to recover the solution xw for wells and applying
// xw to update Well State
void recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const;
int numWells() const;
/// return true if wells are available in the reservoir
bool wellsActive() const;
void setWellsActive(const bool wells_active);
/// return true if wells are available on this process
bool localWellsActive() const;
bool getWellConvergence(Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const;
/// upate the dynamic lists related to economic limits
void updateListEconLimited(const Schedule& schedule,
const int current_step,
const Wells* wells_struct,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
WellCollection* wellCollection() const;
protected:
bool wells_active_;
const Wells* wells_;
const std::vector< const Well* > wells_ecl_;
// the number of wells in this process
// trying not to use things from Wells struct
// TODO: maybe a better name to emphasize it is local?
const int number_of_wells_;
const int number_of_phases_;
using WellInterfacePtr = std::unique_ptr<WellInterface<TypeTag> >;
// a vector of all the wells.
// eventually, the wells_ above should be gone.
// the name is just temporary
// later, might make share_ptr const later.
std::vector<WellInterfacePtr > well_container_;
using ConvergenceReport = typename WellInterface<TypeTag>::ConvergenceReport;
// create the well container
static std::vector<WellInterfacePtr > createWellContainer(const Wells* wells,
const std::vector<const Well*>& wells_ecl,
const int time_step);
// Well collection is used to enforce the group control
WellCollection* well_collection_;
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;
int current_timeIdx_;
PhaseUsage phase_usage_;
std::vector<bool> active_;
const RateConverterType& rate_converter_;
// the number of the cells in the local grid
int number_of_cells_;
long int global_nc_;
// used to better efficiency of calcuation
mutable BVector scaleAddRes_;
void updateWellControls(WellState& xw) const;
void updateGroupControls(WellState& well_state) const;
// setting the well_solutions_ based on well_state.
void updatePrimaryVariables(const WellState& well_state) const;
void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
void computeRepRadiusPerfLength(const Grid& grid);
void computeAverageFormationFactor(Simulator& ebosSimulator,
std::vector<double>& B_avg) const;
void applyVREPGroupControl(WellState& well_state) const;
void computeWellVoidageRates(const WellState& well_state,
std::vector<double>& well_voidage_rates,
std::vector<double>& voidage_conversion_coeffs) const;
// Calculating well potentials for each well
// TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
const std::vector<double>& wellPerfEfficiencyFactors() const;
void calculateEfficiencyFactors();
void computeWellConnectionPressures(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.
int numComponents() const
{
if (numPhases() == 2) {
@@ -141,286 +252,21 @@ enum WellVariablePositions {
return numComp;
}
int numPhases() const;
SimulatorReport assemble(Simulator& ebosSimulator,
const int iterationIdx,
const double dt,
WellState& well_state);
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
void resetWellControlFromState(const WellState& xw) const;
void assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells);
bool only_wells) const;
void
getMobility(const Simulator& ebosSimulator,
const int w,
const int perf,
const int cell_idx,
std::vector<EvalWell>& mob) const;
bool allow_cross_flow(const int w, const Simulator& ebosSimulator) const;
void localInvert(Mat& istlA) const;
void print(Mat& istlA) const;
// substract Binv(D)rw from r;
void apply( BVector& r) const;
// subtract B*inv(D)*C * x from A*x
void apply(const BVector& x, BVector& Ax) const;
// apply well model with scaling of alpha
void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
// xw = inv(D)*(rw - C*x)
void recoverVariable(const BVector& x, BVector& xw) const;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowToEbosPvIdx( const int flowPv ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
std::vector<double>
extractPerfData(const std::vector<double>& in) const;
int numPhases() const;
int numCells() const;
void resetWellControlFromState(const WellState& xw) const;
const Wells& wells() const;
const Wells* wellsPointer() const;
/// return true if wells are available in the reservoir
bool wellsActive() const;
void setWellsActive(const bool wells_active);
/// return true if wells are available on this process
bool localWellsActive() const;
int numWellVars() const;
/// Density of each well perforation
const std::vector<double>& wellPerforationDensities() const;
/// Diff to bhp for each well perforation.
const std::vector<double>& wellPerforationPressureDiffs() const;
EvalWell extendEval(const Eval& in) const;
void setWellVariables(const WellState& xw);
void print(const EvalWell& in) const;
void computeAccumWells();
void computeWellFlux(const int& w, const double& Tw, const IntensiveQuantities& intQuants, const std::vector<EvalWell>& mob_perfcells_dense,
const EvalWell& bhp, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const;
SimulatorReport solveWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state);
void printIf(const int c, const double x, const double y, const double eps, const std::string type) const;
std::vector<double> residual() const;
bool getWellConvergence(Simulator& ebosSimulator,
const int iteration) const;
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw);
void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf) const;
void updateWellState(const BVector& dwells,
WellState& well_state) const;
void updateWellControls(WellState& xw) const;
/// upate the dynamic lists related to economic limits
void updateListEconLimited(const Schedule& schedule,
const int current_step,
const Wells* wells_struct,
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
void computeWellConnectionDensitesPressures(const WellState& xw,
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,
const std::vector<double>& depth_perf,
const double grav);
// Calculating well potentials for each well
// TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials) const;
// TODO: some preparation work, mostly related to group control and RESV,
// 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);
WellCollection* wellCollection() const;
const std::vector<double>&
wellPerfEfficiencyFactors() const;
void calculateEfficiencyFactors();
void computeWellVoidageRates(const WellState& well_state,
std::vector<double>& well_voidage_rates,
std::vector<double>& voidage_conversion_coeffs) const;
void applyVREPGroupControl(WellState& well_state) const;
protected:
bool wells_active_;
const Wells* wells_;
const std::vector< const Well* > wells_ecl_;
// Well collection is used to enforce the group control
WellCollection* well_collection_;
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;
int current_timeIdx_;
PhaseUsage phase_usage_;
std::vector<bool> active_;
const VFPProperties* vfp_properties_;
double gravity_;
const RateConverterType& rate_converter_;
// The efficiency factor for each connection. It is specified based on wells and groups,
// We calculate the factor for each connection for the computation of contributions to the mass balance equations.
// By default, they should all be one.
std::vector<double> well_perforation_efficiency_factors_;
// the depth of the all the cell centers
// for standard Wells, it the same with the perforation depth
std::vector<double> cell_depths_;
std::vector<double> pv_;
std::vector<double> well_perforation_densities_;
std::vector<double> well_perforation_pressure_diffs_;
std::vector<double> wpolymer_;
std::vector<double> wsolvent_;
std::vector<double> wells_rep_radius_;
std::vector<double> wells_perf_length_;
std::vector<double> wells_bore_diameter_;
std::vector<EvalWell> wellVariables_;
std::vector<double> F0_;
Mat duneB_;
Mat duneC_;
Mat invDuneD_;
BVector resWell_;
long int global_nc_;
mutable BVector Cx_;
mutable BVector invDrw_;
mutable BVector scaleAddRes_;
double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
double dWellFractionMax() const {return param_.dwell_fraction_max_; }
// protected methods
EvalWell getBhp(const int wellIdx) const;
EvalWell getQs(const int wellIdx, const int compIdx) const;
EvalWell wellVolumeFraction(const int wellIdx, const int compIdx) const;
EvalWell wellVolumeFractionScaled(const int wellIdx, const int compIdx) const;
// Q_p / (Q_w + Q_g + Q_o) for three phase cases.
EvalWell wellSurfaceVolumeFraction(const int well_index, const int compIdx) const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
const int well_number) const;
using WellMapType = typename WellState::WellMapType;
using WellMapEntryType = typename WellState::mapentry_t;
// a tuple type for ratio limit check.
// first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three
// values should not be used.
// second value indicates whehter there is only one connection left.
// third value indicates the indx of the worst-offending connection.
// the last value indicates the extent of the violation for the worst-offending connection, which is defined by
// the ratio of the actual value to the value of the violated limit.
using RatioCheckTuple = std::tuple<bool, bool, int, double>;
enum ConnectionIndex {
INVALIDCONNECTION = -10000
};
RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
const WellMapEntryType& map_entry) const;
RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
const WellMapEntryType& map_entry) const;
void updateWellStateWithTarget(const WellControls* wc,
const int current,
const int well_index,
WellState& xw) const;
bool wellHasTHPConstraints(const int well_index) const;
// TODO: maybe we should provide a light version of computeWellFlux, which does not include the
// calculation of the derivatives
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const int well_index,
std::vector<double>& well_flux) const;
double mostStrictBhpFromBhpLimits(const int well_index) const;
// TODO: maybe it should be improved to be calculate general rates for THP control later
std::vector<double>
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const int well_index,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential) const;
double wsolvent(const int well_index) const;
double wpolymer(const int well_index) const;
void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const;
void computeRepRadiusPerfLength(const Grid& grid);
};

2881
opm/autodiff/StandardWellsDense_impl.hpp
View File

File diff suppressed because it is too large Load Diff

11
opm/autodiff/WellHelpers.hpp
View File

@@ -35,7 +35,7 @@ namespace Opm {
namespace wellhelpers
{
inline
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
const int well,
@@ -147,6 +147,15 @@ namespace Opm {
return dp;
}
inline
double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
const double rho, const double gravity) {
const double dh = vfp_ref_depth - well_ref_depth;
const double dp = rho * gravity * dh;
return dp;
}
template <class Vector>
inline
Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,

313
opm/autodiff/WellInterface.hpp Normal file
View File

@@ -0,0 +1,313 @@
/*
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_WELLINTERFACE_HEADER_INCLUDED
#define OPM_WELLINTERFACE_HEADER_INCLUDED
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/VFPInjProperties.hpp>
#include <opm/autodiff/VFPProdProperties.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
#include<dune/common/fmatrix.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/matrixmatrix.hh>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <string>
#include <memory>
#include <vector>
#include <cassert>
namespace Opm
{
template<typename TypeTag>
class WellInterface
{
public:
using WellState = WellStateFullyImplicitBlackoilDense;
typedef BlackoilModelParameters ModelParameters;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
static const int numEq = BlackoilIndices::numEq;
typedef double Scalar;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
/// Constructor
WellInterface(const Well* well, const int time_step, const Wells* wells);
/// Virutal destructor
virtual ~WellInterface() {}
/// Well name.
const std::string& name() const;
/// Well type, INJECTOR or PRODUCER.
WellType wellType() const;
/// Well controls
WellControls* wellControls() const;
void setVFPProperties(const VFPProperties* vfp_properties_arg);
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 = 0;
/// a struct to collect information about the convergence checking
struct ConvergenceReport {
struct ProblemWell {
std::string well_name;
std::string phase_name;
};
bool converged = true;
bool nan_residual_found = false;
std::vector<ProblemWell> nan_residual_wells;
// We consider Inf is large residual here
bool too_large_residual_found = false;
std::vector<ProblemWell> too_large_residual_wells;
ConvergenceReport& operator+=(const ConvergenceReport& rhs) {
converged = converged && rhs.converged;
nan_residual_found = nan_residual_found || rhs.nan_residual_found;
if (rhs.nan_residual_found) {
for (const ProblemWell& well : rhs.nan_residual_wells) {
nan_residual_wells.push_back(well);
}
}
too_large_residual_found = too_large_residual_found || rhs.too_large_residual_found;
if (rhs.too_large_residual_found) {
for (const ProblemWell& well : rhs.too_large_residual_wells) {
too_large_residual_wells.push_back(well);
}
}
return *this;
}
};
virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const = 0;
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state) = 0;
virtual void assembleWellEq(Simulator& ebosSimulator,
const double dt,
WellState& well_state,
bool only_wells) = 0;
void updateListEconLimited(const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
void setWellEfficiencyFactor(const double efficiency_factor);
void computeRepRadiusPerfLength(const Grid& grid, const std::map<int, int>& cartesian_to_compressed);
/// 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,
WellState& well_state) const = 0;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const = 0;
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const = 0;
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;
virtual void updateWellStateWithTarget(const int current,
WellState& xw) const = 0;
virtual void updateWellControl(WellState& xw,
wellhelpers::WellSwitchingLogger& logger) const = 0;
virtual void updatePrimaryVariables(const WellState& well_state) const = 0;
protected:
// to indicate a invalid connection
static const int INVALIDCONNECTION = -100000;
const Well* well_ecl_;
const int current_step_;
// the index of well in Wells struct
int index_of_well_;
// well type
// INJECTOR or PRODUCER
enum WellType well_type_;
// number of phases
int number_of_phases_;
// component fractions for each well
// typically, it should apply to injection wells
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_;
// reference depth for the BHP
double ref_depth_;
double well_efficiency_factor_;
// cell index for each well perforation
std::vector<int> well_cells_;
// saturation table nubmer for each well perforation
std::vector<int> saturation_table_number_;
// representative radius of the perforations, used in shear calculation
std::vector<double> perf_rep_radius_;
// length of the perforations, use in shear calculation
std::vector<double> perf_length_;
// well bore diameter
std::vector<double> bore_diameters_;
const PhaseUsage* phase_usage_;
bool getAllowCrossFlow() const;
const std::vector<bool>* active_;
const VFPProperties* vfp_properties_;
double gravity_;
const std::vector<bool>& active() const;
const PhaseUsage& phaseUsage() const;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowToEbosPvIdx( const int flowPv ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
// TODO: it is dumplicated with StandardWellsDense
int numComponents() const;
double wsolvent() const;
double wpolymer() const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const;
bool wellHasTHPConstraints() const;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
double mostStrictBhpFromBhpLimits() const;
// a tuple type for ratio limit check.
// first value indicates whether ratio limit is violated, when the ratio limit is not violated, the following three
// values should not be used.
// second value indicates whehter there is only one connection left.
// third value indicates the indx of the worst-offending connection.
// the last value indicates the extent of the violation for the worst-offending connection, which is defined by
// the ratio of the actual value to the value of the violated limit.
using RatioCheckTuple = std::tuple<bool, bool, int, double>;
RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const;
RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const;
};
}
#include "WellInterface_impl.hpp"
#endif // OPM_WELLINTERFACE_HEADER_INCLUDED

768
opm/autodiff/WellInterface_impl.hpp Normal file
View File

@@ -0,0 +1,768 @@
/*
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/>.
*/
namespace Opm
{
template<typename TypeTag>
WellInterface<TypeTag>::
WellInterface(const Well* well, const int time_step, const Wells* wells)
: well_ecl_(well)
, current_step_(time_step)
{
if (!well) {
OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface");
}
if (time_step < 0) {
OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
}
if (!wells) {
OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface");
}
const std::string& well_name = well->name();
// looking for the location of the well in the wells struct
int index_well;
for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
if (well_name == std::string(wells->name[index_well])) {
break;
}
}
// should not enter the constructor if the well does not exist in the wells struct
// here, just another assertion.
assert(index_well != wells->number_of_wells);
index_of_well_ = index_well;
well_type_ = wells->type[index_well];
number_of_phases_ = wells->number_of_phases;
// copying the comp_frac
{
comp_frac_.resize(number_of_phases_);
const int index_begin = index_well * number_of_phases_;
std::copy(wells->comp_frac + index_begin,
wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
}
well_controls_ = wells->ctrls[index_well];
ref_depth_ = wells->depth_ref[index_well];
// perforations related
{
const int perf_index_begin = wells->well_connpos[index_well];
const int perf_index_end = wells->well_connpos[index_well + 1];
number_of_perforations_ = perf_index_end - perf_index_begin;
first_perf_ = perf_index_begin;
well_cells_.resize(number_of_perforations_);
std::copy(wells->well_cells + perf_index_begin,
wells->well_cells + perf_index_end,
well_cells_.begin() );
well_index_.resize(number_of_perforations_);
std::copy(wells->WI + perf_index_begin,
wells->WI + perf_index_end,
well_index_.begin() );
saturation_table_number_.resize(number_of_perforations_);
std::copy(wells->sat_table_id + perf_index_begin,
wells->sat_table_id + perf_index_end,
saturation_table_number_.begin() );
}
well_efficiency_factor_ = 1.0;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
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 */)
{
phase_usage_ = phase_usage_arg;
active_ = active_arg;
gravity_ = gravity_arg;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
setVFPProperties(const VFPProperties* vfp_properties_arg)
{
vfp_properties_ = vfp_properties_arg;
}
template<typename TypeTag>
const std::string&
WellInterface<TypeTag>::
name() const
{
return well_ecl_->name();
}
template<typename TypeTag>
WellType
WellInterface<TypeTag>::
wellType() const
{
return well_type_;
}
template<typename TypeTag>
WellControls*
WellInterface<TypeTag>::
wellControls() const
{
return well_controls_;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
getAllowCrossFlow() const
{
return well_ecl_->getAllowCrossFlow();
}
template<typename TypeTag>
const std::vector<bool>&
WellInterface<TypeTag>::
active() const
{
assert(active_);
return *active_;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
setWellEfficiencyFactor(const double efficiency_factor)
{
well_efficiency_factor_ = efficiency_factor;
}
template<typename TypeTag>
const PhaseUsage&
WellInterface<TypeTag>::
phaseUsage() const
{
assert(phase_usage_);
return *phase_usage_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
if (phaseIdx > 2 )
return phaseIdx;
return phaseToComp[ phaseIdx ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[ 3 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx
};
if (flowPv > 2 )
return flowPv;
return flowToEbos[ flowPv ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
assert(phaseIdx < 3);
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
return flowToEbos[ phaseIdx ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
numComponents() const
{
// TODO: how about two phase polymer
if (number_of_phases_ == 2) {
return 2;
}
int numComp = FluidSystem::numComponents;
if (has_solvent) {
numComp++;
}
return numComp;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
wsolvent() const
{
if (!has_solvent) {
return 0.0;
}
WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
if (injection.injectorType == WellInjector::GAS) {
double solvent_fraction = well_ecl_->getSolventFraction(current_step_);
return solvent_fraction;
}
assert(false);
return 0.0;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
wpolymer() const
{
if (!has_polymer) {
return 0.0;
}
WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_);
if (injection.injectorType == WellInjector::WATER) {
const double polymer_injection_concentration = polymer.m_polymerConcentration;
return polymer_injection_concentration;
}
assert(false); // TODO: find a more logical way to handle this situation
return 0.0;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
mostStrictBhpFromBhpLimits() const
{
double bhp;
// initial bhp value, making the value not usable
switch( well_type_ ) {
case INJECTOR:
bhp = std::numeric_limits<double>::max();
break;
case PRODUCER:
bhp = -std::numeric_limits<double>::max();
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
}
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
// finding a BHP constraint
if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
// get the bhp constraint value, it should always be postive assummingly
const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
switch(well_type_) {
case INJECTOR: // using the lower bhp contraint from Injectors
if (bhp_target < bhp) {
bhp = bhp_target;
}
break;
case PRODUCER:
if (bhp_target > bhp) {
bhp = bhp_target;
}
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
} // end of switch
}
}
return bhp;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
wellHasTHPConstraints() const
{
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
return true;
}
}
return false;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
const Opm::PhaseUsage& pu = *phase_usage_;
const int np = number_of_phases_;
if (econ_production_limits.onMinOilRate()) {
assert(active()[Oil]);
const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
const double min_oil_rate = econ_production_limits.minOilRate();
if (std::abs(oil_rate) < min_oil_rate) {
return true;
}
}
if (econ_production_limits.onMinGasRate() ) {
assert(active()[Gas]);
const double gas_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
const double min_gas_rate = econ_production_limits.minGasRate();
if (std::abs(gas_rate) < min_gas_rate) {
return true;
}
}
if (econ_production_limits.onMinLiquidRate() ) {
assert(active()[Oil]);
assert(active()[Water]);
const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate;
const double min_liquid_rate = econ_production_limits.minLiquidRate();
if (std::abs(liquid_rate) < min_liquid_rate) {
return true;
}
}
if (econ_production_limits.onMinReservoirFluidRate()) {
OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
}
return false;
}
template<typename TypeTag>
typename WellInterface<TypeTag>::RatioCheckTuple
WellInterface<TypeTag>::
checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
bool water_cut_limit_violated = false;
int worst_offending_connection = INVALIDCONNECTION;
bool last_connection = false;
double violation_extent = -1.0;
const int np = number_of_phases_;
const Opm::PhaseUsage& pu = *phase_usage_;
const int well_number = index_of_well_;
assert(active()[Oil]);
assert(active()[Water]);
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate;
double water_cut;
if (std::abs(liquid_rate) != 0.) {
water_cut = water_rate / liquid_rate;
} else {
water_cut = 0.0;
}
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
if (water_cut > max_water_cut_limit) {
water_cut_limit_violated = true;
}
if (water_cut_limit_violated) {
// need to handle the worst_offending_connection
const int perf_start = first_perf_;
const int perf_number = number_of_perforations_;
std::vector<double> water_cut_perf(perf_number);
for (int perf = 0; perf < perf_number; ++perf) {
const int i_perf = perf_start + perf;
const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
if (std::abs(liquid_perf_rate) != 0.) {
water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
} else {
water_cut_perf[perf] = 0.;
}
}
last_connection = (perf_number == 1);
if (last_connection) {
worst_offending_connection = 0;
violation_extent = water_cut_perf[0] / max_water_cut_limit;
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
double max_water_cut_perf = 0.;
for (int perf = 0; perf < perf_number; ++perf) {
if (water_cut_perf[perf] > max_water_cut_perf) {
worst_offending_connection = perf;
max_water_cut_perf = water_cut_perf[perf];
}
}
assert(max_water_cut_perf != 0.);
assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
violation_extent = max_water_cut_perf / max_water_cut_limit;
}
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
template<typename TypeTag>
typename WellInterface<TypeTag>::RatioCheckTuple
WellInterface<TypeTag>::
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
// TODO: not sure how to define the worst-offending connection when more than one
// ratio related limit is violated.
// The defintion used here is that we define the violation extent based on the
// ratio between the value and the corresponding limit.
// For each violated limit, we decide the worst-offending connection separately.
// Among the worst-offending connections, we use the one has the biggest violation
// extent.
bool any_limit_violated = false;
bool last_connection = false;
int worst_offending_connection = INVALIDCONNECTION;
double violation_extent = -1.0;
if (econ_production_limits.onMaxWaterCut()) {
const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state);
bool water_cut_violated = std::get<0>(water_cut_return);
if (water_cut_violated) {
any_limit_violated = true;
const double violation_extent_water_cut = std::get<3>(water_cut_return);
if (violation_extent_water_cut > violation_extent) {
violation_extent = violation_extent_water_cut;
worst_offending_connection = std::get<2>(water_cut_return);
last_connection = std::get<1>(water_cut_return);
}
}
}
if (econ_production_limits.onMaxGasOilRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
}
if (econ_production_limits.onMaxWaterGasRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
}
if (econ_production_limits.onMaxGasLiquidRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
}
if (any_limit_violated) {
assert(worst_offending_connection >=0);
assert(violation_extent > 1.);
}
return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
updateListEconLimited(const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const
{
// economic limits only apply for production wells.
if (wellType() != PRODUCER) {
return;
}
// flag to check if the mim oil/gas rate limit is violated
bool rate_limit_violated = false;
const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_);
// if no limit is effective here, then continue to the next well
if ( !econ_production_limits.onAnyEffectiveLimit() ) {
return;
}
const std::string well_name = name();
// for the moment, we only handle rate limits, not handling potential limits
// the potential limits should not be difficult to add
const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
if (quantity_limit == WellEcon::POTN) {
const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n")
+ std::string("All the limits will be evaluated based on RATE. ");
OpmLog::warning("NOT_SUPPORTING_POTN", msg);
}
if (econ_production_limits.onAnyRateLimit()) {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state);
}
if (rate_limit_violated) {
if (econ_production_limits.endRun()) {
const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n")
+ std::string("the program will keep running after ") + well_name + std::string(" is closed");
OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message);
}
if (econ_production_limits.validFollowonWell()) {
OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
}
if (well_ecl_->getAutomaticShutIn()) {
list_econ_limited.addShutWell(well_name);
const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit");
OpmLog::info(msg);
} else {
list_econ_limited.addStoppedWell(well_name);
const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit");
OpmLog::info(msg);
}
// the well is closed, not need to check other limits
return;
}
// checking for ratio related limits, mostly all kinds of ratio.
bool ratio_limits_violated = false;
RatioCheckTuple ratio_check_return;
if (econ_production_limits.onAnyRatioLimit()) {
ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state);
ratio_limits_violated = std::get<0>(ratio_check_return);
}
if (ratio_limits_violated) {
const bool last_connection = std::get<1>(ratio_check_return);
const int worst_offending_connection = std::get<2>(ratio_check_return);
assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_));
const int cell_worst_offending_connection = well_cells_[worst_offending_connection];
list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ")
+ well_name + std::string(" will be closed due to economic limit");
OpmLog::info(msg);
if (last_connection) {
list_econ_limited.addShutWell(well_name);
const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed");
OpmLog::info(msg2);
}
}
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
computeRepRadiusPerfLength(const Grid& grid,
const std::map<int, int>& cartesian_to_compressed)
{
const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
const int nperf = number_of_perforations_;
perf_rep_radius_.clear();
perf_length_.clear();
bore_diameters_.clear();
perf_rep_radius_.resize(nperf);
perf_length_.resize(nperf);
bore_diameters_.resize(nperf);
// COMPDAT handling
const auto& completionSet = well_ecl_->getCompletions(current_step_);
for (size_t c=0; c<completionSet.size(); c++) {
const auto& completion = completionSet.get(c);
if (completion.getState() == WellCompletion::OPEN) {
const int i = completion.getI();
const int j = completion.getJ();
const int k = completion.getK();
const int* cpgdim = cart_dims;
const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
const std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << name() << ')');
}
const int cell = cgit->second;
{
double radius = 0.5*completion.getDiameter();
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
WellCompletion::DirectionEnum direction = completion.getDirection();
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
break;
case Opm::WellCompletion::DirectionEnum::Y:
re = std::sqrt(cubical[0] * cubical[2] / M_PI);
perf_length = cubical[1];
break;
case Opm::WellCompletion::DirectionEnum::Z:
re = std::sqrt(cubical[0] * cubical[1] / M_PI);
perf_length = cubical[2];
break;
default:
OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
}
const double repR = std::sqrt(re * radius);
perf_rep_radius_.push_back(repR);
perf_length_.push_back(perf_length);
bore_diameters_.push_back(2. * radius);
}
}
}
}
}

101
opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
View File

@@ -108,101 +108,6 @@ namespace Opm
}
}
}
// TODO: the reason to keep this is to avoid getting defaulted value BHP
// some facilities needed from opm-parser or opm-core
// It is a little tricky, since sometimes before applying group control, the only
// available constraints in the well_controls is the defaulted BHP value, and it
// is really not desirable to use this value to enter the Newton iterations.
setWellSolutions(pu);
}
/// Set wellSolutions() based on the base class members.
void setWellSolutions(const PhaseUsage& pu)
{
// Set nw and np, or return if no wells.
if (wells_.get() == nullptr) {
return;
}
const int nw = wells_->number_of_wells;
if (nw == 0) {
return;
}
const int np = wells_->number_of_phases;
const int numComp = pu.has_solvent? np+1:np;
well_solutions_.clear();
well_solutions_.resize(nw * numComp, 0.0);
std::vector<double> g = {1.0,1.0,0.01};
for (int w = 0; w < nw; ++w) {
WellControls* wc = wells_->ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = currentControls()[w];
well_controls_set_current( wc, current);
const WellType& well_type = wells_->type[w];
switch (well_controls_iget_type(wc, current)) {
case THP: // Intentional fall-through
case BHP:
if (well_type == INJECTOR) {
for (int p = 0; p < np; ++p) {
well_solutions_[w] += wellRates()[np*w + p] * wells_->comp_frac[np*w + p];
}
} else {
for (int p = 0; p < np; ++p) {
well_solutions_[w] += g[p] * wellRates()[np*w + p];
}
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
wellSolutions()[w] = bhp()[w];
break;
}
double total_rates = 0.0;
for (int p = 0; p < np; ++p) {
total_rates += g[p] * wellRates()[np*w + p];
}
const int waterpos = pu.phase_pos[Water];
const int gaspos = pu.phase_pos[Gas];
assert(np > 2 || (np == 2 && !pu.phase_used[Gas]));
// assumes the gas fractions are stored after water fractions
if(std::abs(total_rates) > 0) {
if( pu.phase_used[Water] ) {
wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates;
}
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ;
}
if( pu.has_solvent) {
wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates;
}
} else {
if( pu.phase_used[Water] ) {
wellSolutions()[nw + w] = wells_->comp_frac[np*w + waterpos];
}
if( pu.phase_used[Gas] ) {
wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos];
}
if (pu.has_solvent) {
wellSolutions()[3*nw + w] = 0;
}
}
}
}
@@ -212,11 +117,6 @@ namespace Opm
init(wells, state.pressure(), dummy_state, pu) ;
}
/// One rate per phase and well connection.
std::vector<double>& wellSolutions() { return well_solutions_; }
const std::vector<double>& wellSolutions() const { return well_solutions_; }
/// One rate pr well connection.
std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
@@ -252,7 +152,6 @@ namespace Opm
private:
std::vector<double> well_solutions_;
std::vector<double> perfRateSolvent_;
};

58
tests/TESTWELLMODEL.DATA Normal file
View File

@@ -0,0 +1,58 @@
RUNSPEC
DIMENS
5 5 4 /
OIL
GAS
WATER
GRID
DX
100*100. /
DY
100*50. /
DZ
100*10. /
TOPS
25*2500 /
PORO
100*0.3 /
PERMX
100*10. /
PERMY
100*20. /
PERMZ
100*1. /
SCHEDULE
WELSPECS
'PROD1' 'P' 5 5 2525 'OIL' /
'INJE1' 'I' 1 1 2505 'GAS' /
/
COMPDAT
'PROD1' 5 5 3 4 'OPEN' 2* 0.15 /
'INJE1' 1 1 1 4 'OPEN' 2* 0.15 /
/
WCONINJE
'INJE1' 'WATER' 'OPEN' 'RATE' 1000. /
/
WCONPROD
'PROD1' 'OPEN' 'GRAT' 2* 50000. /
/
TSTEP
10 /
END

11
tests/test_multisegmentwells.cpp
View File

@@ -83,15 +83,6 @@ struct SetupMSW {
std::unique_ptr<GridInit> grid_init(new GridInit(ecl_state, porv));
const Grid& grid = grid_init->grid();
// Create material law manager.
std::vector<int> compressed_to_cartesianIdx;
Opm::createGlobalCellArray(grid, compressed_to_cartesianIdx);
std::shared_ptr<MaterialLawManager> material_law_manager(new MaterialLawManager());
material_law_manager->initFromDeck(deck, ecl_state, compressed_to_cartesianIdx);
std::unique_ptr<FluidProps> fluidprops(new FluidProps(deck, ecl_state, material_law_manager, grid));
const size_t current_timestep = 0;
// dummy_dynamic_list_econ_lmited
@@ -115,7 +106,7 @@ struct SetupMSW {
std::unordered_set<std::string>());
const Wells* wells = wells_manager.c_wells();
const auto wells_ecl = ecl_state.getSchedule().getWells(current_timestep);
const auto& wells_ecl = ecl_state.getSchedule().getWells(current_timestep);
ms_wells.reset(new Opm::MultisegmentWells(wells, &(wells_manager.wellCollection()), wells_ecl, current_timestep));
};

194
tests/test_wellmodel.cpp Normal file
View File

@@ -0,0 +1,194 @@
/*
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/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE WellModelTest
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/filesystem.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/DynamicListEconLimited.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <ebos/eclproblem.hh>
#include <ewoms/common/start.hh>
#include <opm/autodiff/StandardWell.hpp>
#include <opm/autodiff/StandardWellsDense.hpp>
// maybe should just include BlackoilModelEbos.hpp
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
}
}
using StandardWell = Opm::StandardWell<TTAG(EclFlowProblem)>;
struct SetupTest {
using Grid = UnstructuredGrid;
using GridInit = Opm::GridInit<Grid>;
SetupTest ()
{
Opm::ParseContext parse_context;
Opm::Parser parser;
auto deck = parser.parseFile("TESTWELLMODEL.DATA", parse_context);
ecl_state = std::make_unique<const Opm::EclipseState>(deck , parse_context);
// Create grid.
const std::vector<double>& porv =
ecl_state->get3DProperties().getDoubleGridProperty("PORV").getData();
std::unique_ptr<GridInit> grid_init(new GridInit(*ecl_state, porv));
const Grid& grid = grid_init->grid();
// Create material law manager.
std::vector<int> compressed_to_cartesianIdx;
Opm::createGlobalCellArray(grid, compressed_to_cartesianIdx);
// dummy_dynamic_list_econ_lmited
const Opm::DynamicListEconLimited dummy_dynamic_list;
current_timestep = 0;
// Create wells.
wells_manager = std::make_unique<const Opm::WellsManager> (*ecl_state,
current_timestep,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::dimensions(grid),
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
dummy_dynamic_list,
false,
std::unordered_set<std::string>());
};
std::unique_ptr<const Opm::WellsManager> wells_manager;
std::unique_ptr<const Opm::EclipseState> ecl_state;
int current_timestep;
};
BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
SetupTest setup_test;
const Wells* wells = setup_test.wells_manager->c_wells();
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);
}
BOOST_AUTO_TEST_CASE(TestBehavoir) {
SetupTest setup_test;
const Wells* wells_struct = setup_test.wells_manager->c_wells();
const auto& wells_ecl = setup_test.ecl_state->getSchedule().getWells(setup_test.current_timestep);
const int current_timestep = setup_test.current_timestep;
std::vector<std::unique_ptr<const StandardWell> > wells;
{
const int nw = wells_struct ? (wells_struct->number_of_wells) : 0;
for (int w = 0; w < nw; ++w) {
const std::string well_name(wells_struct->name[w]);
size_t index_well = 0;
for (; index_well < wells_ecl.size(); ++index_well) {
if (well_name == wells_ecl[index_well]->name()) {
break;
}
}
// 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) );
}
}
// first well, it is a production well from the deck
{
const auto& well = wells[0];
BOOST_CHECK_EQUAL(well->name(), "PROD1");
BOOST_CHECK(well->wellType() == PRODUCER);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numWellEq == 3);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);
const auto& control = well_controls_get_current(wc);
BOOST_CHECK(control >= 0);
// GAS RATE CONTROL
const auto& distr = well_controls_iget_distr(wc, control);
BOOST_CHECK(distr[0] == 0.);
BOOST_CHECK(distr[1] == 0.);
BOOST_CHECK(distr[2] == 1.);
}
// second well, it is the injection well from the deck
{
const auto& well = wells[1];
BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->wellType() == INJECTOR);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numWellEq == 3);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);
const auto& control = well_controls_get_current(wc);
BOOST_CHECK(control >= 0);
// WATER RATE CONTROL
const auto& distr = well_controls_iget_distr(wc, control);
BOOST_CHECK(distr[0] == 1.);
BOOST_CHECK(distr[1] == 0.);
BOOST_CHECK(distr[2] == 0.);
}
}