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/parser/eclipse/EclipseState/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 #include <dune/istl/umfpack.hh>
36 
37 #include <array>
38 #include <memory>
39 
40 namespace Opm
41 {
42 
43 class ConvergenceReport;
44 class GroupState;
45 class Schedule;
46 class WellContributions;
47 template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
48 class WellState;
49 
50 template<typename FluidSystem, typename Indices, typename Scalar>
52 {
53 public:
54 #if HAVE_CUDA || HAVE_OPENCL
56  void addWellContribution(WellContributions& wellContribs) const;
57 #endif
58 
59 protected:
60  // TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
61 
62  // TODO: we need to have order for the primary variables and also the order for the well equations.
63  // sometimes, they are similar, while sometimes, they can have very different forms.
64 
65  // Table showing the primary variable indices, depending on what phases are present:
66  //
67  // WOG OG WG WO W/O/G (single phase)
68  // GTotal 0 0 0 0 0
69  // WFrac 1 -1000 1 1 -1000
70  // GFrac 2 1 -1000 -1000 -1000
71  // Spres 3 2 2 2 1
72 
73  static constexpr bool has_water = (Indices::waterSaturationIdx >= 0);
74  static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
75  static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 0;
76 
77  // In the implementation, one should use has_wfrac_variable
78  // rather than has_water to check if you should do something
79  // with the variable at the WFrac location, similar for GFrac.
80  static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1;
81  static constexpr bool has_gfrac_variable = has_gas && has_oil;
82 
83  static constexpr int GTotal = 0;
84  static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
85  static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
86  static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1;
87 
88  // the number of well equations TODO: it should have a more general strategy for it
89  static constexpr int numWellEq = Indices::numPhases + 1;
90 
91  // sparsity pattern for the matrices
92  // [A C^T [x = [ res
93  // B D ] x_well] res_well]
94 
95  // the vector type for the res_well and x_well
96  using VectorBlockWellType = Dune::FieldVector<Scalar, numWellEq>;
97  using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
98 
99  using VectorBlockType = Dune::FieldVector<Scalar, Indices::numEq>;
100  using BVector = Dune::BlockVector<VectorBlockType>;
101 
102  // the matrix type for the diagonal matrix D
103  using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, numWellEq>;
104  using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
105 
106  // the matrix type for the non-diagonal matrix B and C^T
107  using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, Indices::numEq>;
108  using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
109 
110  // TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
111  // EvalR (Eval), EvalW, EvalRW
112  // TODO: for now, we only use one type to save some implementation efforts, while improve later.
113  using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
114  using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
115 
117 
118  void initMatrixAndVectors(const int num_cells) const;
119  void initPrimaryVariablesEvaluation() const;
120 
121  void assembleControlEq(const WellState& well_state,
122  const GroupState& group_state,
123  const Schedule& schedule,
124  const SummaryState& summaryState,
125  const Well::InjectionControls& inj_controls,
126  const Well::ProductionControls& prod_controls,
127  const double rho,
128  DeferredLogger& deferred_logger);
129 
130  void assembleDefaultPressureEq(const int seg,
131  WellState& well_state) const;
132 
133 
134  // assemble pressure equation for ICD segments
135  void assembleICDPressureEq(const int seg,
136  const UnitSystem& unit_system,
137  WellState& well_state,
138  DeferredLogger& deferred_logger) const;
139 
140 
141  void assemblePressureEq(const int seg,
142  const UnitSystem& unit_system,
143  WellState& well_state,
144  DeferredLogger& deferred_logger) const;
145 
146  void checkConvergenceControlEq(const WellState& well_state,
147  ConvergenceReport& report,
148  const double tolerance_pressure_ms_wells,
149  const double tolerance_wells,
150  const double max_residual_allowed,
151  DeferredLogger& deferred_logger) const;
152 
155  const std::vector<double>& B_avg,
156  DeferredLogger& deferred_logger,
157  const double max_residual_allowed,
158  const double tolerance_wells,
159  const double relaxed_inner_tolerance_flow_ms_well,
160  const double tolerance_pressure_ms_wells,
161  const double relaxed_inner_tolerance_pressure_ms_well,
162  const bool relax_tolerance) const;
163 
164  // handling the overshooting and undershooting of the fractions
165  void processFractions(const int seg) const;
166 
167  // xw = inv(D)*(rw - C*x)
168  void recoverSolutionWell(const BVector& x,
169  BVectorWell& xw) const;
170 
171  void updatePrimaryVariables(const WellState& well_state) const;
172 
173  void updateUpwindingSegments();
174 
175  // updating the well_state based on well solution dwells
176  void updateWellState(const BVectorWell& dwells,
177  const double relaxation_factor,
178  const double DFLimit,
179  const double max_pressure_change) const;
180 
181  void computeSegmentFluidProperties(const EvalWell& temperature,
182  const EvalWell& saltConcentration,
183  int pvt_region_index);
184 
185  EvalWell getBhp() const;
186  EvalWell getFrictionPressureLoss(const int seg) const;
187  EvalWell getHydroPressureLoss(const int seg) const;
188  EvalWell getQs(const int comp_idx) const;
189  EvalWell getSegmentGTotal(const int seg) const;
190  EvalWell getSegmentPressure(const int seg) const;
191  EvalWell getSegmentRate(const int seg, const int comp_idx) const;
192  EvalWell getSegmentRateUpwinding(const int seg,
193  const size_t comp_idx) const;
194  EvalWell getSegmentSurfaceVolume(const EvalWell& temperature,
195  const EvalWell& saltConcentration,
196  const int pvt_region_index,
197  const int seg_idx) const;
198  EvalWell getWQTotal() const;
199 
200 
201  std::vector<Scalar> getWellResiduals(const std::vector<Scalar>& B_avg,
202  DeferredLogger& deferred_logger) const;
203 
204  double getControlTolerance(const WellState& well_state,
205  const double tolerance_wells,
206  const double tolerance_pressure_ms_wells,
207  DeferredLogger& deferred_logger) const;
208 
209  double getResidualMeasureValue(const WellState& well_state,
210  const std::vector<double>& residuals,
211  const double tolerance_wells,
212  const double tolerance_pressure_ms_wells,
213  DeferredLogger& deferred_logger) const;
214 
215  void handleAccelerationPressureLoss(const int seg,
216  WellState& well_state) const;
217 
218  // pressure drop for Autonomous ICD segment (WSEGAICD)
219  EvalWell pressureDropAutoICD(const int seg,
220  const UnitSystem& unit_system) const;
221 
222  // pressure drop for Spiral ICD segment (WSEGSICD)
223  EvalWell pressureDropSpiralICD(const int seg) const;
224 
225  // pressure drop for sub-critical valve (WSEGVALV)
226  EvalWell pressureDropValve(const int seg) const;
227 
228  void updateThp(WellState& well_state,
229  const double rho,
230  DeferredLogger& deferred_logger) const;
231 
232  void updateWellStateFromPrimaryVariables(WellState& well_state,
233  const double rho,
234  DeferredLogger& deferred_logger) const;
235 
236  // fraction value of the primary variables
237  // should we just use member variables to store them instead of calculating them again and again
238  EvalWell volumeFraction(const int seg,
239  const unsigned compIdx) const;
240 
241  // F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
242  EvalWell volumeFractionScaled(const int seg,
243  const int comp_idx) const;
244 
245  // basically Q_p / \sigma_p Q_p
246  EvalWell surfaceVolumeFraction(const int seg,
247  const int comp_idx) const;
248 
249  // convert a Eval from reservoir to contain the derivative related to wells
250  EvalWell extendEval(const Eval& in) const;
251 
253 
254  // TODO, the following should go to a class for computing purpose
255  // two off-diagonal matrices
256  mutable OffDiagMatWell duneB_;
257  mutable OffDiagMatWell duneC_;
258  // "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets.
259  mutable DiagMatWell duneD_;
260 
264  mutable std::shared_ptr<Dune::UMFPack<DiagMatWell> > duneDSolver_;
265 
266  // residuals of the well equations
267  mutable BVectorWell resWell_;
268 
269  // the values for the primary varibles
270  // based on different solutioin strategies, the wells can have different primary variables
271  mutable std::vector<std::array<double, numWellEq> > primary_variables_;
272 
273  // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
274  mutable std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
275 
276  // the upwinding segment for each segment based on the flow direction
277  std::vector<int> upwinding_segments_;
278 
279  // the densities of segment fluids
280  // we should not have this member variable
281  std::vector<EvalWell> segment_densities_;
282 
283  // the mass rate of the segments
284  std::vector<EvalWell> segment_mass_rates_;
285 
286  // the viscosity of the segments
287  std::vector<EvalWell> segment_viscosities_;
288 
289  std::vector<std::vector<EvalWell>> segment_phase_densities_;
290  std::vector<std::vector<EvalWell>> segment_phase_fractions_;
291  std::vector<std::vector<EvalWell>> segment_phase_viscosities_;
292 
293  // depth difference between perforations and the perforated grid cells
294  std::vector<double> cell_perforation_depth_diffs_;
295  // pressure correction due to the different depth of the perforation and
296  // center depth of the grid block
297  std::vector<double> cell_perforation_pressure_diffs_;
298 };
299 
300 }
301 
302 #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:52
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:229
std::shared_ptr< Dune::UMFPack< DiagMatWell > > duneDSolver_
solver for diagonal matrix
Definition: MultisegmentWellEval.hpp:264
Definition: MultisegmentWellGeneric.hpp:41
This class serves to eliminate the need to include the WellContributions into the matrix (with –matri...
Definition: WellContributions.hpp:61
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:26