2017-08-17 08:49:54 -05:00
/*
Copyright 2017 SINTEF Digital , Mathematics and Cybernetics .
Copyright 2017 Statoil ASA .
This file is part of the Open Porous Media project ( OPM ) .
OPM is free software : you can redistribute it and / or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation , either version 3 of the License , or
( at your option ) any later version .
OPM is distributed in the hope that it will be useful ,
but WITHOUT ANY WARRANTY ; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE . See the
GNU General Public License for more details .
You should have received a copy of the GNU General Public License
along with OPM . If not , see < http : //www.gnu.org/licenses/>.
*/
# ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
# define OPM_MULTISEGMENTWELL_HEADER_INCLUDED
# include <opm/autodiff/WellInterface.hpp>
2017-09-13 10:10:34 -05:00
# include <dune/istl/solvers.hh>
2017-08-17 08:49:54 -05:00
namespace Opm
{
template < typename TypeTag >
class MultisegmentWell : public WellInterface < TypeTag >
{
public :
typedef WellInterface < TypeTag > Base ;
// TODO: the WellState does not have any information related to segments
using typename Base : : WellState ;
2017-09-01 10:37:51 -05:00
using typename Base : : Simulator ;
2017-09-08 04:09:54 -05:00
using typename Base : : IntensiveQuantities ;
using typename Base : : FluidSystem ;
2017-09-01 10:37:51 -05:00
using typename Base : : ModelParameters ;
2017-09-08 07:58:57 -05:00
using typename Base : : MaterialLaw ;
2017-08-17 08:49:54 -05:00
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
// TODO: should I begin with the old primary variable or the new fraction based variable systems?
// Let us begin with the new one
2017-09-12 06:16:57 -05:00
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
2017-08-17 08:49:54 -05:00
enum WellVariablePositions {
GTotal = 0 ,
WFrac = 1 ,
GFrac = 2 ,
SPres = 3
} ;
/// the number of well equations // TODO: it should have a more general strategy for it
static const int numWellEq = 4 ;
using typename Base : : Scalar ;
using typename Base : : ConvergenceReport ;
/// the number of reservior equations
using Base : : numEq ;
/// the matrix and vector types for the reservoir
using typename Base : : Mat ;
using typename Base : : BVector ;
using typename Base : : Eval ;
// sparsity pattern for the matrices
// [A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
typedef Dune : : FieldVector < Scalar , numWellEq > VectorBlockWellType ;
typedef Dune : : BlockVector < VectorBlockWellType > BVectorWell ;
// the matrix type for the diagonal matrix D
typedef Dune : : FieldMatrix < Scalar , numWellEq , numWellEq > DiagMatrixBlockWellType ;
typedef Dune : : BCRSMatrix < DiagMatrixBlockWellType > DiagMatWell ;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune : : FieldMatrix < Scalar , numWellEq , numEq > OffDiagMatrixBlockWellType ;
typedef Dune : : BCRSMatrix < OffDiagMatrixBlockWellType > OffDiagMatWell ;
// TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
// EvalR (Eval), EvalW, EvalRW
// TODO: for now, we only use one type to save some implementation efforts, while improve later.
typedef DenseAd : : Evaluation < double , /*size=*/ numEq + numWellEq > EvalWell ;
MultisegmentWell ( const Well * well , const int time_step , const Wells * wells ) ;
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 ,
2017-09-04 10:58:37 -05:00
WellState & well_state ) const ;
2017-08-17 08:49:54 -05:00
/// 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 ;
/// 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 ,
2017-09-05 04:09:58 -05:00
std : : vector < double > & well_potentials ) ;
2017-08-17 08:49:54 -05:00
virtual void updatePrimaryVariables ( const WellState & well_state ) const ;
2017-08-29 10:26:36 -05:00
virtual void solveEqAndUpdateWellState ( const ModelParameters & param ,
WellState & well_state ) ; // const?
2017-09-05 04:20:43 -05:00
virtual void calculateExplicitQuantities ( const Simulator & ebosSimulator ,
const WellState & well_state ) ; // should be const?
2017-09-05 04:09:58 -05:00
2017-08-28 10:30:22 -05:00
/// number of segments for this well
/// int number_of_segments_;
int numberOfSegments ( ) const ;
int numberOfPerforations ( ) const ;
2017-08-17 08:49:54 -05:00
protected :
int number_segments_ ;
// components of the pressure drop to be included
WellSegment : : CompPressureDropEnum compPressureDrop ( ) const ;
// multi-phase flow model
WellSegment : : MultiPhaseModelEnum multiphaseModel ( ) const ;
// get the SegmentSet from the well_ecl_
const SegmentSet & segmentSet ( ) const ;
2017-09-04 05:42:23 -05:00
// protected member variables from the Base class
2017-08-17 08:49:54 -05:00
using Base : : well_ecl_ ;
using Base : : number_of_perforations_ ; // TODO: can use well_ecl_?
2017-09-01 10:37:51 -05:00
using Base : : current_step_ ;
2017-09-04 05:42:23 -05:00
using Base : : index_of_well_ ;
using Base : : number_of_phases_ ;
2017-08-17 08:49:54 -05:00
using Base : : well_cells_ ; // TODO: are the perforation orders same with StandardWell or Wells?
using Base : : well_index_ ;
2017-09-04 05:42:23 -05:00
using Base : : well_type_ ;
2017-09-04 07:47:27 -05:00
using Base : : first_perf_ ;
2017-09-08 07:58:57 -05:00
using Base : : saturation_table_number_ ;
using Base : : well_efficiency_factor_ ;
2017-09-12 07:26:32 -05:00
using Base : : gravity_ ;
2017-08-17 08:49:54 -05:00
using Base : : well_controls_ ;
2017-09-04 05:42:23 -05:00
// protected functions from the Base class
using Base : : active ;
2017-09-04 07:47:27 -05:00
using Base : : phaseUsage ;
2017-09-04 10:58:37 -05:00
using Base : : name ;
2017-09-05 08:16:04 -05:00
using Base : : numComponents ;
2017-09-08 04:09:54 -05:00
using Base : : flowToEbosPvIdx ;
using Base : : flowPhaseToEbosPhaseIdx ;
2017-09-08 07:58:57 -05:00
using Base : : flowPhaseToEbosCompIdx ;
using Base : : getAllowCrossFlow ;
2017-09-04 05:42:23 -05:00
2017-08-17 08:49:54 -05:00
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// indices of the gird blocks that segments locate at.
// TODO: the grid cell related to a segment should be calculated based on the location
// of the segment node.
// As the current temporary solution, the grid cell related to a segment determined by the
// first perforation cell related to the segment.
// when no perforation is related to the segment, use it outlet segment's cell.
2017-09-01 10:37:51 -05:00
// TODO: it can be a source of error
2017-08-17 08:49:54 -05:00
std : : vector < int > segment_cell_ ;
2017-08-31 11:22:48 -05:00
// the completions that is related to each segment
// the completions's ids are their location in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
2017-09-01 10:37:51 -05:00
std : : vector < std : : vector < int > > segment_perforations_ ;
// the inlet segments for each segment. It is for convinience and efficiency reason
// the original segment structure is defined as a gathering tree structure based on outlet_segment
2017-09-04 05:42:23 -05:00
// the reason that we can not use the old way of WellOps, which is based on the Eigen matrix and vector.
// TODO: can we use DUNE FieldMatrix and FieldVector.
2017-09-01 10:37:51 -05:00
std : : vector < std : : vector < int > > segment_inlets_ ;
2017-08-31 11:22:48 -05:00
2017-08-17 08:49:54 -05:00
// Things are easy to get from SegmentSet
// segment_volume_, segment_cross_area_, segment_length_(total length), segment_depth_
// segment_internal_diameter_, segment_roughness_
// outlet_segment_., in the outlet_segment, we store the ID of the segment, we will need to use numberToLocation to get
// their location in the segmentSet
// segment number is an ID of the segment, it is specified in the deck
// get the loation of the segment with a segment number in the segmentSet
int numberToLocation ( const int segment_number ) const ;
// TODO, the following should go to a class for computing purpose
// two off-diagonal matrices
mutable OffDiagMatWell duneB_ ;
mutable OffDiagMatWell duneC_ ;
// diagonal matrix for the well
2017-09-12 10:32:27 -05:00
// TODO: if we decided not to invert it, we better change the name of it
2017-09-13 10:48:08 -05:00
mutable DiagMatWell duneD_ ;
2017-08-17 08:49:54 -05:00
// several vector used in the matrix calculation
mutable BVectorWell Bx_ ;
mutable BVector scaleAddRes_ ;
// residuals of the well equations
2017-09-01 10:37:51 -05:00
mutable BVectorWell resWell_ ;
2017-08-17 08:49:54 -05:00
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
// TODO: should we introduce a data structure for segment to simplify this?
2017-08-29 10:26:36 -05:00
// or std::vector<std::vector<double> >
mutable std : : vector < std : : array < double , numWellEq > > primary_variables_ ;
2017-08-17 08:49:54 -05:00
// the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
2017-08-29 10:26:36 -05:00
mutable std : : vector < std : : array < EvalWell , numWellEq > > primary_variables_evaluation_ ;
2017-09-05 08:16:04 -05:00
// pressure correction due to the different depth of the perforation and
// center depth of the grid block
std : : vector < double > perforation_cell_pressure_diffs_ ;
2017-09-08 04:09:54 -05:00
// depth difference between the segment and the peforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std : : vector < double > segment_perforation_depth_diffs_ ;
2017-09-05 08:16:04 -05:00
// the intial component compistion of segments
2017-09-05 08:34:31 -05:00
std : : vector < std : : vector < double > > segment_comp_initial_ ;
2017-09-05 08:16:04 -05:00
2017-09-08 04:09:54 -05:00
// the densities of segment fluids
// TODO: if it turned out it is only used to calculate the pressure difference,
// we should not have this member variable
std : : vector < EvalWell > segment_densities_ ;
2017-09-12 07:26:32 -05:00
std : : vector < double > segment_depth_diffs_ ;
2017-09-01 10:37:51 -05:00
void initMatrixAndVectors ( const int num_cells ) const ;
2017-08-29 10:26:36 -05:00
// protected functions
// EvalWell getBhp(); this one should be something similar to getSegmentPressure();
// EvalWell getQs(); this one should be something similar to getSegmentRates()
// EValWell wellVolumeFractionScaled, wellVolumeFraction, wellSurfaceVolumeFraction ... these should have different names, and probably will be needed.
// bool crossFlowAllowed(const Simulator& ebosSimulator) const; probably will be needed
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell ( const BVector & x , BVectorWell & xw ) const ;
// updating the well_state based on well solution dwells
void updateWellState ( const BVectorWell & dwells ,
const BlackoilModelParameters & param ,
WellState & well_state ) const ;
2017-09-05 04:20:43 -05:00
// computing the accumulation term for later use in well mass equations
2017-09-05 08:16:04 -05:00
void computeInitialComposition ( ) ;
// compute the pressure difference between the perforation and cell center
void computePerfCellPressDiffs ( const Simulator & ebosSimulator ) ;
// fraction value of the primary variables
// should we just use member variables to store them instead of calculating them again and again
EvalWell volumeFraction ( const int seg , const int comp_idx ) const ;
// F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
EvalWell volumeFractionScaled ( const int seg , const int comp_idx ) const ;
// basically Q_p / \sigma_p Q_p
EvalWell surfaceVolumeFraction ( const int seg , const int comp_idx ) const ;
2017-09-05 04:20:43 -05:00
2017-09-08 04:09:54 -05:00
void computePerfRate ( const IntensiveQuantities & int_quants ,
const std : : vector < EvalWell > & mob_perfcells ,
const int seg ,
const double well_index ,
const EvalWell & segment_pressure ,
const bool & allow_cf ,
std : : vector < EvalWell > & cq_s ) const ;
// convert a Eval from reservoir to contain the derivative related to wells
EvalWell extendEval ( const Eval & in ) const ;
// compute the densities of the mixture in the segments
// TODO: probably other fluid properties also.
// They will be treated implicitly, so they need to be of Evaluation type
void computeSegmentFluidProperties ( const Simulator & ebosSimulator ,
const WellState & well_state ) ;
EvalWell getSegmentPressure ( const int seg ) const ;
2017-09-08 07:58:57 -05:00
EvalWell getSegmentRate ( const int seg , const int comp_idx ) const ;
2017-09-11 10:31:42 -05:00
EvalWell getSegmentGTotal ( const int seg ) const ;
2017-09-08 07:58:57 -05:00
// get the mobility for specific perforation
void getMobility ( const Simulator & ebosSimulator ,
const int perf ,
std : : vector < EvalWell > & mob ) const ;
2017-09-11 10:31:42 -05:00
2017-09-15 06:42:15 -05:00
void assembleControlEq ( ) const ;
2017-09-12 06:16:57 -05:00
2017-09-15 06:42:15 -05:00
void assemblePressureEq ( const int seg ) const ;
2017-09-12 07:26:32 -05:00
// hytrostatic pressure loss
EvalWell getHydorPressureLoss ( const int seg ) const ;
2017-09-12 10:13:30 -05:00
// handling the overshooting and undershooting of the fractions
2017-09-12 10:32:27 -05:00
void processFractions ( const int seg ) const ;
void updateWellStateFromPrimaryVariabls ( WellState & well_state ) const ;
2017-08-17 08:49:54 -05:00
} ;
2017-09-13 10:10:34 -05:00
// obtain y = D^-1 * x
template < typename MatrixType , typename VectorType >
VectorType
invDX ( MatrixType D , VectorType x )
{
// TODO: checking the problem related to use reference parameter
VectorType y ( x . size ( ) ) ;
y = 0. ;
Dune : : MatrixAdapter < MatrixType , VectorType , VectorType > linearOperator ( D ) ;
// Sequential incomplete LU decomposition as the preconditioner
Dune : : SeqILU0 < MatrixType , VectorType , VectorType > preconditioner ( D , 1.0 ) ;
2017-09-15 06:42:15 -05:00
// Preconditioned BICGSTAB solver
2017-09-13 10:10:34 -05:00
Dune : : BiCGSTABSolver < VectorType > linsolver ( linearOperator ,
preconditioner ,
1.e-6 , // desired residual reduction factor
50 , // maximum number of iterations
0 ) ; // verbosity of the solver
// Object storing some statistics about the solving process
Dune : : InverseOperatorResult statistics ;
// Solve
linsolver . apply ( y , x , statistics ) ;
return y ;
}
2017-08-17 08:49:54 -05:00
}
# include "MultisegmentWell_impl.hpp"
# endif // OPM_MULTISEGMENTWELL_HEADER_INCLUDED