All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilMultiSegmentModel_impl.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED
21 #define OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED
22 
23 #include <opm/autodiff/BlackoilMultiSegmentModel.hpp>
24 
25 #include <opm/autodiff/AutoDiffBlock.hpp>
26 #include <opm/autodiff/AutoDiffHelpers.hpp>
27 #include <opm/autodiff/GridHelpers.hpp>
28 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
29 #include <opm/autodiff/GeoProps.hpp>
30 #include <opm/autodiff/WellDensitySegmented.hpp>
31 #include <opm/autodiff/VFPProperties.hpp>
32 #include <opm/autodiff/VFPProdProperties.hpp>
33 #include <opm/autodiff/VFPInjProperties.hpp>
34 
35 #include <opm/core/grid.h>
36 #include <opm/core/linalg/LinearSolverInterface.hpp>
37 #include <opm/core/linalg/ParallelIstlInformation.hpp>
38 #include <opm/core/props/rock/RockCompressibility.hpp>
39 #include <opm/common/ErrorMacros.hpp>
40 #include <opm/common/Exceptions.hpp>
41 #include <opm/parser/eclipse/Units/Units.hpp>
42 #include <opm/core/well_controls.h>
43 #include <opm/core/utility/parameters/ParameterGroup.hpp>
44 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
45 
46 #include <cassert>
47 #include <cmath>
48 #include <iostream>
49 #include <iomanip>
50 #include <limits>
51 #include <vector>
52 //#include <fstream>
53 
54 
55 namespace Opm {
56 
57 
58  template <class Grid>
60  BlackoilMultiSegmentModel(const typename Base::ModelParameters& param,
61  const Grid& grid ,
62  const BlackoilPropsAdFromDeck& fluid,
63  const DerivedGeology& geo ,
64  const RockCompressibility* rock_comp_props,
65  const MultisegmentWells& well_model,
66  const NewtonIterationBlackoilInterface& linsolver,
67  std::shared_ptr< const EclipseState > eclState,
68  const bool has_disgas,
69  const bool has_vapoil,
70  const bool terminal_output)
71  : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
72  eclState, has_disgas, has_vapoil, terminal_output)
73  {
74  }
75 
76 
77 
78 
79 
80  template <class Grid>
81  void
84  const ReservoirState& reservoir_state,
85  const WellState& well_state)
86  {
87  const double dt = timer.currentStepLength();
88  pvdt_ = geo_.poreVolume() / dt;
89  if (active_[Gas]) {
90  updatePrimalVariableFromState(reservoir_state);
91  }
92 
93  const int nw = wellsMultiSegment().size();
94 
95  if ( !msWellOps().has_multisegment_wells ) {
96  wellModel().segVDt() = V::Zero(nw);
97  return;
98  }
99 
100  const int nseg_total = well_state.numSegments();
101  std::vector<double> segment_volume;
102  segment_volume.reserve(nseg_total);
103  for (int w = 0; w < nw; ++w) {
104  WellMultiSegmentConstPtr well = wellsMultiSegment()[w];
105  const std::vector<double>& segment_volume_well = well->segmentVolume();
106  segment_volume.insert(segment_volume.end(), segment_volume_well.begin(), segment_volume_well.end());
107  }
108  assert(int(segment_volume.size()) == nseg_total);
109  wellModel().segVDt() = Eigen::Map<V>(segment_volume.data(), nseg_total) / dt;
110  }
111 
112 
113 
114 
115 
116  template <class Grid>
117  void
118  BlackoilMultiSegmentModel<Grid>::makeConstantState(SolutionState& state) const
119  {
120  Base::makeConstantState(state);
121  state.segp = ADB::constant(state.segp.value());
122  state.segqs = ADB::constant(state.segqs.value());
123  }
124 
125 
126 
127 
128 
129  template <class Grid>
130  SimulatorReport
132  assemble(const ReservoirState& reservoir_state,
133  WellState& well_state,
134  const bool initial_assembly)
135  {
136  using namespace Opm::AutoDiffGrid;
137 
138  // TODO: include VFP effect.
139  // If we have VFP tables, we need the well connection
140  // pressures for the "simple" hydrostatic correction
141  // between well depth and vfp table depth.
142  // if (isVFPActive()) {
143  // SolutionState state = asImpl().variableState(reservoir_state, well_state);
144  // SolutionState state0 = state;
145  // asImpl().makeConstantState(state0);
146  // asImpl().computeWellConnectionPressures(state0, well_state);
147  // }
148 
149  // Possibly switch well controls and updating well state to
150  // get reasonable initial conditions for the wells
151  wellModel().updateWellControls(well_state);
152 
153  // TODO: I do not think the multi_segment well can handle group control yet
154  if (asImpl().wellModel().wellCollection()->groupControlActive()) {
155  // enforce VREP control when necessary.
156  Base::applyVREPGroupControl(reservoir_state, well_state);
157 
158  asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
159  }
160 
161  // Create the primary variables.
162  SolutionState state = asImpl().variableState(reservoir_state, well_state);
163 
164  if (initial_assembly) {
165  // Create the (constant, derivativeless) initial state.
166  SolutionState state0 = state;
167  asImpl().makeConstantState(state0);
168  // Compute initial accumulation contributions
169  // and well connection pressures.
170  asImpl().computeAccum(state0, 0);
171  wellModel().computeSegmentFluidProperties(state0);
172  const int np = numPhases();
173  assert(np == int(wellModel().segmentCompSurfVolumeInitial().size()));
174  for (int phase = 0; phase < np; ++phase) {
175  wellModel().segmentCompSurfVolumeInitial()[phase] = wellModel().segmentCompSurfVolumeCurrent()[phase].value();
176  }
177 
178  asImpl().computeWellConnectionPressures(state0, well_state);
179  }
180 
181  // OPM_AD_DISKVAL(state.pressure);
182  // OPM_AD_DISKVAL(state.saturation[0]);
183  // OPM_AD_DISKVAL(state.saturation[1]);
184  // OPM_AD_DISKVAL(state.saturation[2]);
185  // OPM_AD_DISKVAL(state.rs);
186  // OPM_AD_DISKVAL(state.rv);
187  // OPM_AD_DISKVAL(state.qs);
188  // OPM_AD_DISKVAL(state.bhp);
189 
190  // -------- Mass balance equations --------
191  asImpl().assembleMassBalanceEq(state);
192 
193  // -------- Well equations ----------
194  if ( ! wellsActive() ) {
195  SimulatorReport report;
196  return report;
197  }
198 
199  wellModel().computeSegmentFluidProperties(state);
200 
201  const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
202  wellModel().computeSegmentPressuresDelta(gravity);
203 
204  std::vector<ADB> mob_perfcells;
205  std::vector<ADB> b_perfcells;
206  SimulatorReport report;
207  wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
208  if (param_.solve_welleq_initially_ && initial_assembly) {
209  // solve the well equations as a pre-processing step
210  report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
211  }
212 
213  // the perforation flux here are different
214  // it is related to the segment location
215  V aliveWells;
216  std::vector<ADB> cq_s;
217  wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
218  wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
219  wellModel().addWellFluxEq(cq_s, state, residual_);
220  asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
221  wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
222  return report;
223  }
224 
225 
226 
227 
228 
229  template <class Grid>
230  SimulatorReport
231  BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
232  const std::vector<ADB>& b_perfcells,
233  const ReservoirState& reservoir_state,
234  SolutionState& state,
235  WellState& well_state)
236  {
237  SimulatorReport report = Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
238 
239  if (report.converged) {
240  // We must now update the state.segp and state.segqs members,
241  // that the base version does not know about.
242  const int np = numPhases();
243  const int nseg_total =well_state.numSegments();
244  {
245  // We will set the segp primary variable to the new ones,
246  // but we do not change the derivatives here.
247  ADB::V new_segp = Eigen::Map<ADB::V>(well_state.segPress().data(), nseg_total);
248  // Avoiding the copy below would require a value setter method
249  // in AutoDiffBlock.
250  std::vector<ADB::M> old_segp_derivs = state.segp.derivative();
251  state.segp = ADB::function(std::move(new_segp), std::move(old_segp_derivs));
252  }
253  {
254  // Need to reshuffle well rates, from phase running fastest
255  // to wells running fastest.
256  // The transpose() below switches the ordering.
257  const DataBlock segrates = Eigen::Map<const DataBlock>(well_state.segPhaseRates().data(), nseg_total, np).transpose();
258  ADB::V new_segqs = Eigen::Map<const V>(segrates.data(), nseg_total * np);
259  std::vector<ADB::M> old_segqs_derivs = state.segqs.derivative();
260  state.segqs = ADB::function(std::move(new_segqs), std::move(old_segqs_derivs));
261  }
262 
263  // This is also called by the base version, but since we have updated
264  // state.segp we must call it again.
265  asImpl().computeWellConnectionPressures(state, well_state);
266  }
267 
268  return report;
269  }
270 
271 
272 
273 
274  template <class Grid>
275  void
276  BlackoilMultiSegmentModel<Grid>::
277  computeWellConnectionPressures(const SolutionState& state,
278  const WellState& well_state)
279  {
280  const int np = numPhases();
281  const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
282  std::vector<ADB> fluid_density(np, ADB::null());
283  // TODO: make sure the order of the density and the order of the kr are the same.
284  for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
285  const int canonicalPhaseIdx = canph_[phaseIdx];
286  fluid_density[phaseIdx] = fluidDensity(canonicalPhaseIdx, sd_.rq[phaseIdx].b, state.rs, state.rv);
287  }
288  wellModel().computeWellConnectionPressures(state, well_state, kr_adb, fluid_density);
289  }
290 
291 } // namespace Opm
292 
293 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
BlackoilMultiSegmentModel(const typename Base::ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const MultisegmentWells &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output)
Construct the model.
Definition: BlackoilMultiSegmentModel_impl.hpp:60
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
Definition: BlackoilMultiSegmentModel.hpp:35
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilMultiSegmentModel_impl.hpp:83
virtual double currentStepLength() const =0
Current step length.
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilMultiSegmentModel_impl.hpp:132
A model implementation for three-phase black oil with support for multi-segment wells.
Definition: BlackoilMultiSegmentModel.hpp:55
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
Class for handling the multi-segment well model.
Definition: MultisegmentWells.hpp:55