My Project
MultisegmentWellEval.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED
23 #define OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED
24 
25 #include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
26 
27 #include <opm/material/densead/Evaluation.hpp>
28 
29 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
30 
31 #include <dune/common/fmatrix.hh>
32 #include <dune/common/fvector.hh>
33 #include <dune/istl/bcrsmatrix.hh>
34 #include <dune/istl/bvector.hh>
35 
36 #include <array>
37 #include <memory>
38 
39 namespace Dune {
40 template<class Matrix> class UMFPack;
41 }
42 
43 namespace Opm
44 {
45 
46 class ConvergenceReport;
47 class GroupState;
48 class Schedule;
49 class WellContributions;
50 template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
51 class WellState;
52 
53 template<typename FluidSystem, typename Indices, typename Scalar>
55 {
56 public:
58  void addWellContribution(WellContributions& wellContribs) const;
59 
60 protected:
61  // TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
62 
63  // TODO: we need to have order for the primary variables and also the order for the well equations.
64  // sometimes, they are similar, while sometimes, they can have very different forms.
65 
66  // Table showing the primary variable indices, depending on what phases are present:
67  //
68  // WOG OG WG WO W/O/G (single phase)
69  // WQTotal 0 0 0 0 0
70  // WFrac 1 -1000 1 1 -1000
71  // GFrac 2 1 -1000 -1000 -1000
72  // Spres 3 2 2 2 1
73 
74  static constexpr bool has_water = (Indices::waterSaturationIdx >= 0);
75  static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
76  static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0;
77 
78  // In the implementation, one should use has_wfrac_variable
79  // rather than has_water to check if you should do something
80  // with the variable at the WFrac location, similar for GFrac.
81  static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1;
82  static constexpr bool has_gfrac_variable = has_gas && has_oil;
83 
84  static constexpr int WQTotal = 0;
85  static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
86  static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
87  static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1;
88 
89  // the number of well equations TODO: it should have a more general strategy for it
90  static constexpr int numWellEq = Indices::numPhases + 1;
91 
92  // sparsity pattern for the matrices
93  // [A C^T [x = [ res
94  // B D ] x_well] res_well]
95 
96  // the vector type for the res_well and x_well
97  using VectorBlockWellType = Dune::FieldVector<Scalar, numWellEq>;
98  using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
99 
100  using VectorBlockType = Dune::FieldVector<Scalar, Indices::numEq>;
101  using BVector = Dune::BlockVector<VectorBlockType>;
102 
103  // the matrix type for the diagonal matrix D
104  using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, numWellEq>;
105  using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
106 
107  // the matrix type for the non-diagonal matrix B and C^T
108  using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, Indices::numEq>;
109  using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
110 
111  // TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
112  // EvalR (Eval), EvalW, EvalRW
113  // TODO: for now, we only use one type to save some implementation efforts, while improve later.
114  using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
115  using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
116 
118 
119  void initMatrixAndVectors(const int num_cells) const;
120  void initPrimaryVariablesEvaluation() const;
121 
122  void assembleControlEq(const WellState& well_state,
123  const GroupState& group_state,
124  const Schedule& schedule,
125  const SummaryState& summaryState,
126  const Well::InjectionControls& inj_controls,
127  const Well::ProductionControls& prod_controls,
128  const double rho,
129  DeferredLogger& deferred_logger);
130 
131  void assembleDefaultPressureEq(const int seg,
132  WellState& well_state) const;
133 
134 
135  // assemble pressure equation for ICD segments
136  void assembleICDPressureEq(const int seg,
137  const UnitSystem& unit_system,
138  WellState& well_state,
139  DeferredLogger& deferred_logger) const;
140 
141 
142  void assemblePressureEq(const int seg,
143  const UnitSystem& unit_system,
144  WellState& well_state,
145  DeferredLogger& deferred_logger) const;
146 
147  void checkConvergenceControlEq(const WellState& well_state,
148  ConvergenceReport& report,
149  const double tolerance_pressure_ms_wells,
150  const double tolerance_wells,
151  const double max_residual_allowed,
152  DeferredLogger& deferred_logger) const;
153 
156  const std::vector<double>& B_avg,
157  DeferredLogger& deferred_logger,
158  const double max_residual_allowed,
159  const double tolerance_wells,
160  const double relaxed_inner_tolerance_flow_ms_well,
161  const double tolerance_pressure_ms_wells,
162  const double relaxed_inner_tolerance_pressure_ms_well,
163  const bool relax_tolerance) const;
164 
165  // handling the overshooting and undershooting of the fractions
166  void processFractions(const int seg) const;
167 
168  // xw = inv(D)*(rw - C*x)
169  void recoverSolutionWell(const BVector& x,
170  BVectorWell& xw) const;
171 
172  void updatePrimaryVariables(const WellState& well_state) const;
173 
174  void updateUpwindingSegments();
175 
176  // updating the well_state based on well solution dwells
177  void updatePrimaryVariablesNewton(const BVectorWell& dwells,
178  const double relaxation_factor,
179  const double DFLimit,
180  const double max_pressure_change) const;
181 
182  void computeSegmentFluidProperties(const EvalWell& temperature,
183  const EvalWell& saltConcentration,
184  int pvt_region_index,
185  DeferredLogger& deferred_logger);
186 
187  EvalWell getBhp() const;
188  EvalWell getFrictionPressureLoss(const int seg) const;
189  EvalWell getHydroPressureLoss(const int seg) const;
190  EvalWell getQs(const int comp_idx) const;
191  EvalWell getSegmentWQTotal(const int seg) const;
192  EvalWell getSegmentPressure(const int seg) const;
193  EvalWell getSegmentRate(const int seg, const int comp_idx) const;
194  EvalWell getSegmentRateUpwinding(const int seg,
195  const size_t comp_idx) const;
196  EvalWell getSegmentSurfaceVolume(const EvalWell& temperature,
197  const EvalWell& saltConcentration,
198  const int pvt_region_index,
199  const int seg_idx) const;
200  EvalWell getWQTotal() const;
201 
202 
203  std::pair<bool, std::vector<Scalar> >
204  getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
205  DeferredLogger& deferred_logger) const;
206 
207  double getControlTolerance(const WellState& well_state,
208  const double tolerance_wells,
209  const double tolerance_pressure_ms_wells,
210  DeferredLogger& deferred_logger) const;
211 
212  double getResidualMeasureValue(const WellState& well_state,
213  const std::vector<double>& residuals,
214  const double tolerance_wells,
215  const double tolerance_pressure_ms_wells,
216  DeferredLogger& deferred_logger) const;
217 
218  void handleAccelerationPressureLoss(const int seg,
219  WellState& well_state) const;
220 
221  // pressure drop for Autonomous ICD segment (WSEGAICD)
222  EvalWell pressureDropAutoICD(const int seg,
223  const UnitSystem& unit_system) const;
224 
225  // pressure drop for Spiral ICD segment (WSEGSICD)
226  EvalWell pressureDropSpiralICD(const int seg) const;
227 
228  // pressure drop for sub-critical valve (WSEGVALV)
229  EvalWell pressureDropValve(const int seg) const;
230 
231  void updateThp(WellState& well_state,
232  const double rho,
233  DeferredLogger& deferred_logger) const;
234 
235  void updateWellStateFromPrimaryVariables(WellState& well_state,
236  const double rho,
237  DeferredLogger& deferred_logger) const;
238 
239  // fraction value of the primary variables
240  // should we just use member variables to store them instead of calculating them again and again
241  EvalWell volumeFraction(const int seg,
242  const unsigned compIdx) const;
243 
244  // F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
245  EvalWell volumeFractionScaled(const int seg,
246  const int comp_idx) const;
247 
248  // basically Q_p / \sigma_p Q_p
249  EvalWell surfaceVolumeFraction(const int seg,
250  const int comp_idx) const;
251 
252  // convert a Eval from reservoir to contain the derivative related to wells
253  EvalWell extendEval(const Eval& in) const;
254 
256 
257  // TODO, the following should go to a class for computing purpose
258  // two off-diagonal matrices
259  mutable OffDiagMatWell duneB_;
260  mutable OffDiagMatWell duneC_;
261  // "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets.
262  mutable DiagMatWell duneD_;
263 
267  mutable std::shared_ptr<Dune::UMFPack<DiagMatWell> > duneDSolver_;
268 
269  // residuals of the well equations
270  mutable BVectorWell resWell_;
271 
272  // the values for the primary varibles
273  // based on different solutioin strategies, the wells can have different primary variables
274  mutable std::vector<std::array<double, numWellEq> > primary_variables_;
275 
276  // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
277  mutable std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
278 
279  // the upwinding segment for each segment based on the flow direction
280  std::vector<int> upwinding_segments_;
281 
282  // the densities of segment fluids
283  // we should not have this member variable
284  std::vector<EvalWell> segment_densities_;
285 
286  // the mass rate of the segments
287  std::vector<EvalWell> segment_mass_rates_;
288 
289  // the viscosity of the segments
290  std::vector<EvalWell> segment_viscosities_;
291 
292  std::vector<std::vector<EvalWell>> segment_phase_densities_;
293  std::vector<std::vector<EvalWell>> segment_phase_fractions_;
294  std::vector<std::vector<EvalWell>> segment_phase_viscosities_;
295 
296  // depth difference between perforations and the perforated grid cells
297  std::vector<double> cell_perforation_depth_diffs_;
298  // pressure correction due to the different depth of the perforation and
299  // center depth of the grid block
300  std::vector<double> cell_perforation_pressure_diffs_;
301 };
302 
303 }
304 
305 #endif // OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition: ConvergenceReport.hpp:36
Definition: DeferredLogger.hpp:57
Definition: GroupState.hpp:34
Definition: MultisegmentWellEval.hpp:55
void addWellContribution(WellContributions &wellContribs) const
add the contribution (C, D, B matrices) of this Well to the WellContributions object
Definition: MultisegmentWellEval.cpp:1827
ConvergenceReport getWellConvergence(const WellState &well_state, const std::vector< double > &B_avg, DeferredLogger &deferred_logger, const double max_residual_allowed, const double tolerance_wells, const double relaxed_inner_tolerance_flow_ms_well, const double tolerance_pressure_ms_wells, const double relaxed_inner_tolerance_pressure_ms_well, const bool relax_tolerance) const
check whether the well equations get converged for this well
Definition: MultisegmentWellEval.cpp:231
std::shared_ptr< Dune::UMFPack< DiagMatWell > > duneDSolver_
solver for diagonal matrix
Definition: MultisegmentWellEval.hpp:267
Definition: MultisegmentWellGeneric.hpp:42
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition: WellContributions.hpp:53
Definition: WellInterfaceIndices.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27