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 #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  DeferredLogger& deferred_logger);
185 
186  EvalWell getBhp() const;
187  EvalWell getFrictionPressureLoss(const int seg) const;
188  EvalWell getHydroPressureLoss(const int seg) const;
189  EvalWell getQs(const int comp_idx) const;
190  EvalWell getSegmentGTotal(const int seg) const;
191  EvalWell getSegmentPressure(const int seg) const;
192  EvalWell getSegmentRate(const int seg, const int comp_idx) const;
193  EvalWell getSegmentRateUpwinding(const int seg,
194  const size_t comp_idx) const;
195  EvalWell getSegmentSurfaceVolume(const EvalWell& temperature,
196  const EvalWell& saltConcentration,
197  const int pvt_region_index,
198  const int seg_idx) const;
199  EvalWell getWQTotal() const;
200 
201 
202  std::pair<bool, std::vector<Scalar> >
203  getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
204  DeferredLogger& deferred_logger) const;
205 
206  double getControlTolerance(const WellState& well_state,
207  const double tolerance_wells,
208  const double tolerance_pressure_ms_wells,
209  DeferredLogger& deferred_logger) const;
210 
211  double getResidualMeasureValue(const WellState& well_state,
212  const std::vector<double>& residuals,
213  const double tolerance_wells,
214  const double tolerance_pressure_ms_wells,
215  DeferredLogger& deferred_logger) const;
216 
217  void handleAccelerationPressureLoss(const int seg,
218  WellState& well_state) const;
219 
220  // pressure drop for Autonomous ICD segment (WSEGAICD)
221  EvalWell pressureDropAutoICD(const int seg,
222  const UnitSystem& unit_system) const;
223 
224  // pressure drop for Spiral ICD segment (WSEGSICD)
225  EvalWell pressureDropSpiralICD(const int seg) const;
226 
227  // pressure drop for sub-critical valve (WSEGVALV)
228  EvalWell pressureDropValve(const int seg) const;
229 
230  void updateThp(WellState& well_state,
231  const double rho,
232  DeferredLogger& deferred_logger) const;
233 
234  void updateWellStateFromPrimaryVariables(WellState& well_state,
235  const double rho,
236  DeferredLogger& deferred_logger) const;
237 
238  // fraction value of the primary variables
239  // should we just use member variables to store them instead of calculating them again and again
240  EvalWell volumeFraction(const int seg,
241  const unsigned compIdx) const;
242 
243  // F_p / g_p, the basic usage of this value is because Q_p = G_t * F_p / G_p
244  EvalWell volumeFractionScaled(const int seg,
245  const int comp_idx) const;
246 
247  // basically Q_p / \sigma_p Q_p
248  EvalWell surfaceVolumeFraction(const int seg,
249  const int comp_idx) const;
250 
251  // convert a Eval from reservoir to contain the derivative related to wells
252  EvalWell extendEval(const Eval& in) const;
253 
255 
256  // TODO, the following should go to a class for computing purpose
257  // two off-diagonal matrices
258  mutable OffDiagMatWell duneB_;
259  mutable OffDiagMatWell duneC_;
260  // "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets.
261  mutable DiagMatWell duneD_;
262 
266  mutable std::shared_ptr<Dune::UMFPack<DiagMatWell> > duneDSolver_;
267 
268  // residuals of the well equations
269  mutable BVectorWell resWell_;
270 
271  // the values for the primary varibles
272  // based on different solutioin strategies, the wells can have different primary variables
273  mutable std::vector<std::array<double, numWellEq> > primary_variables_;
274 
275  // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
276  mutable std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
277 
278  // the upwinding segment for each segment based on the flow direction
279  std::vector<int> upwinding_segments_;
280 
281  // the densities of segment fluids
282  // we should not have this member variable
283  std::vector<EvalWell> segment_densities_;
284 
285  // the mass rate of the segments
286  std::vector<EvalWell> segment_mass_rates_;
287 
288  // the viscosity of the segments
289  std::vector<EvalWell> segment_viscosities_;
290 
291  std::vector<std::vector<EvalWell>> segment_phase_densities_;
292  std::vector<std::vector<EvalWell>> segment_phase_fractions_;
293  std::vector<std::vector<EvalWell>> segment_phase_viscosities_;
294 
295  // depth difference between perforations and the perforated grid cells
296  std::vector<double> cell_perforation_depth_diffs_;
297  // pressure correction due to the different depth of the perforation and
298  // center depth of the grid block
299  std::vector<double> cell_perforation_pressure_diffs_;
300 };
301 
302 }
303 
304 #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:266
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