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
39namespace Dune {
40template<class Matrix> class UMFPack;
41}
42
43namespace Opm
44{
45
46class ConvergenceReport;
47class GroupState;
48class Schedule;
49class WellContributions;
50template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
51class WellState;
52
53template<typename FluidSystem, typename Indices, typename Scalar>
55{
56public:
58 void addWellContribution(WellContributions& wellContribs) const;
59
60protected:
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