BlackoilPressureModel.hpp
1 /*
2  Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
3  Copyright 2016 Statoil AS.
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 #ifndef OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
22 #define OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
23 
24 
25 #include <opm/autodiff/BlackoilModelBase.hpp>
26 #include <opm/core/simulator/BlackoilState.hpp>
27 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
28 #include <opm/autodiff/BlackoilModelParameters.hpp>
29 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
30 
31 #include <algorithm>
32 
33 namespace Opm {
34 
40  template<class Grid, class WellModel>
41  class BlackoilPressureModel : public BlackoilModelBase<Grid, WellModel, BlackoilPressureModel<Grid, WellModel> >
42  {
43  public:
44 
46  friend Base;
47 
48  typedef typename Base::ReservoirState ReservoirState;
49  typedef typename Base::WellState WellState;
50  typedef typename Base::SolutionState SolutionState;
51  typedef typename Base::V V;
52 
67  BlackoilPressureModel(const typename Base::ModelParameters& param,
68  const Grid& grid,
69  const BlackoilPropsAdFromDeck& fluid,
70  const DerivedGeology& geo,
71  const RockCompressibility* rock_comp_props,
72  const StandardWells& std_wells,
73  const NewtonIterationBlackoilInterface& linsolver,
74  std::shared_ptr< const EclipseState> eclState,
75  const bool has_disgas,
76  const bool has_vapoil,
77  const bool terminal_output)
78  : Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
79  eclState, has_disgas, has_vapoil, terminal_output),
80  state0_(3),
81  max_dp_rel_(std::numeric_limits<double>::infinity()),
82  scaling_{ ADB::null(), ADB::null(), ADB::null() }
83  {
84  }
85 
88  const ReservoirState& reservoir_state,
89  const WellState& well_state)
90  {
91  asImpl().wellModel().setStoreWellPerforationFluxesFlag(true);
92  Base::prepareStep(timer, reservoir_state, well_state);
93  max_dp_rel_ = std::numeric_limits<double>::infinity();
94  state0_ = asImpl().variableState(reservoir_state, well_state);
95  asImpl().makeConstantState(state0_);
96  }
97 
98 
102  {
103  // We make a reduced residual object which just
104  // contains the pressure residual.
105  // TODO: directly contruct that object in residual_ instead.
106  const int n1 = residual_.material_balance_eq[0].size();
107  const int n2 = residual_.well_flux_eq.size() + residual_.well_eq.size();
108  const int n_full = residual_.sizeNonLinear();
109  const auto& mb = residual_.material_balance_eq;
110  const auto& fe = residual_.well_flux_eq;
111  const auto& we = residual_.well_eq;
112  LinearisedBlackoilResidual pressure_res = {
113  {
114  // TODO: handle general 2-phase etc.
115  ADB::function(mb[0].value(), { mb[0].derivative()[0], mb[0].derivative()[3], mb[0].derivative()[4] })
116  },
117  ADB::function(fe.value(), { fe.derivative()[0], fe.derivative()[3], fe.derivative()[4] }),
118  ADB::function(we.value(), { we.derivative()[0], we.derivative()[3], we.derivative()[4] }),
119  residual_.matbalscale,
120  residual_.singlePrecision
121  };
122  assert(pressure_res.sizeNonLinear() == n1 + n2);
123  V dx_pressure = linsolver_.computeNewtonIncrement(pressure_res);
124  assert(dx_pressure.size() == n1 + n2);
125  V dx_full = V::Zero(n_full);
126  dx_full.topRows(n1) = dx_pressure.topRows(n1);
127  dx_full.bottomRows(n2) = dx_pressure.bottomRows(n2);
128  return dx_full;
129  }
130 
131  using Base::numPhases;
132  using Base::numMaterials;
133  using Base::wellModel;
134 
135  protected:
136  using Base::asImpl;
137  using Base::linsolver_;
138  using Base::residual_;
139  using Base::sd_;
140  using Base::grid_;
141  using Base::ops_;
142  using Base::has_vapoil_;
143  using Base::has_disgas_;
144 
145  SolutionState state0_;
146  double max_dp_rel_ = std::numeric_limits<double>::infinity();
147  ADB scaling_[3] = { ADB::null(), ADB::null(), ADB::null() };
148 
149 
150 
151 
152  SimulatorReport
153  assemble(const ReservoirState& reservoir_state,
154  WellState& well_state,
155  const bool initial_assembly)
156  {
157  SimulatorReport report;
158 
159  report += Base::assemble(reservoir_state, well_state, initial_assembly);
160  if (initial_assembly) {
161  }
162 
163  // Compute pressure residual.
164  ADB pressure_residual = ADB::constant(V::Zero(residual_.material_balance_eq[0].size()));
165  for (int phase = 0; phase < numPhases(); ++phase) {
166  pressure_residual += residual_.material_balance_eq[phase] * scaling_[phase];
167  }
168  residual_.material_balance_eq[0] = pressure_residual; // HACK
169 
170  // Compute total reservoir volume flux.
171  const int n = sd_.rq[0].mflux.size();
172  V flux = V::Zero(n);
173  for (int phase = 0; phase < numPhases(); ++phase) {
174  UpwindSelector<double> upwind(grid_, ops_, sd_.rq[phase].dh.value());
175  flux += sd_.rq[phase].mflux.value() / upwind.select(sd_.rq[phase].b.value());
176  }
177 
178  // Storing the fluxes in the assemble() method is a bit of
179  // a hack, but alternatives either require a more
180  // significant redesign of the base class or are
181  // inefficient.
182  ReservoirState& s = const_cast<ReservoirState&>(reservoir_state);
183  s.faceflux().resize(n);
184  std::copy_n(flux.data(), n, s.faceflux().begin());
185  if (asImpl().localWellsActive()) {
186  const V& wflux = asImpl().wellModel().getStoredWellPerforationFluxes();
187  assert(int(well_state.perfRates().size()) == wflux.size());
188  std::copy_n(wflux.data(), wflux.size(), well_state.perfRates().begin());
189  }
190 
191  return report;
192  }
193 
194 
195 
196 
197 
198  SolutionState
199  variableState(const ReservoirState& x,
200  const WellState& xw) const
201  {
202  // As Base::variableState(), except making Sw and Xvar constants.
203  std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
204  std::vector<ADB> vars = ADB::variables(vars0);
205  const std::vector<int> indices = asImpl().variableStateIndices();
206  vars[indices[Sw]] = ADB::constant(vars[indices[Sw]].value());
207  vars[indices[Xvar]] = ADB::constant(vars[indices[Xvar]].value());
208  // OpmLog::debug("Injector pressure: " + std::to_string(vars[indices[Bhp]].value()[1]));
209  return asImpl().variableStateExtractVars(x, indices, vars);
210  }
211 
212 
213 
214 
215 
216  void computeAccum(const SolutionState& state,
217  const int aix )
218  {
219  if (aix == 0) {
220  Base::computeAccum(state0_, aix);
221  } else {
222  Base::computeAccum(state, aix);
223  }
224  }
225 
226 
227 
228 
229 
230  void assembleMassBalanceEq(const SolutionState& state)
231  {
232  Base::assembleMassBalanceEq(state);
233 
234  // Compute the scaling factors.
235  // Scaling factors are:
236  // 1/bw, 1/bo - rs/bg, 1/bg - rv/bo
237  assert(numPhases() == 3);
238  assert(numMaterials() == 3);
239  V one = V::Constant(state.pressure.size(), 1.0);
240  scaling_[Water] = one / sd_.rq[Water].b;
241  scaling_[Oil] = one / sd_.rq[Oil].b - state.rs / sd_.rq[Gas].b;
242  scaling_[Gas] = one / sd_.rq[Gas].b - state.rv / sd_.rq[Oil].b;
243  if (has_disgas_ && has_vapoil_) {
244  ADB r_factor = one / (one - state.rs * state.rv);
245  scaling_[Oil] = r_factor * scaling_[Oil];
246  scaling_[Gas] = r_factor * scaling_[Gas];
247  }
248  // @TODO: multiply the oil and gas scale by 1/(1-rs*rv)?
249  // OpmLog::debug("gas scaling: " + std::to_string(scaling_[2].value()[0]));
250  }
251 
252 
253 
254 
255 
256  void updateState(const V& dx,
257  ReservoirState& reservoir_state,
258  WellState& well_state)
259  {
260  // Naively, rs and rv can get overwritten, so we
261  // avoid that by storing.
262  std::vector<double> rs_old = reservoir_state.gasoilratio();
263  std::vector<double> rv_old = reservoir_state.rv();
264  auto hs_old = reservoir_state.hydroCarbonState();
265  auto phasecond_old = Base::phaseCondition_;
266  auto isRs_old = Base::isRs_;
267  auto isRv_old = Base::isRv_;
268  auto isSg_old = Base::isSg_;
269 
270  // Compute the pressure range.
271  const auto minmax_iters = std::minmax_element(reservoir_state.pressure().begin(),
272  reservoir_state.pressure().end());
273  const double range = *minmax_iters.second - *minmax_iters.first;
274 
275  // Use the base class' updateState().
276  Base::updateState(dx, reservoir_state, well_state);
277 
278  // Compute relative change.
279  max_dp_rel_ = dx.head(reservoir_state.pressure().size()).abs().maxCoeff() / range;
280 
281  // Restore rs and rv, also various state flags.
282  reservoir_state.gasoilratio() = rs_old;
283  reservoir_state.rv() = rv_old;
284  reservoir_state.hydroCarbonState() = hs_old;
285  Base::phaseCondition_ = phasecond_old;
286  Base::isRs_ = isRs_old;
287  Base::isRv_ = isRv_old;
288  Base::isSg_ = isSg_old;
289  }
290 
291 
292 
293 
294 
295  bool getConvergence(const SimulatorTimerInterface& /* timer */, const int iteration)
296  {
297  const double tol_p = 1e-10;
298  const double resmax = residual_.material_balance_eq[0].value().abs().maxCoeff();
300  // Only rank 0 does print to std::cout
301  if (iteration == 0) {
302  OpmLog::info("\nIter Res(p) Delta(p)\n");
303  }
304  std::ostringstream os;
305  os.precision(3);
306  os.setf(std::ios::scientific);
307  os << std::setw(4) << iteration;
308  os << std::setw(11) << resmax;
309  os << std::setw(11) << max_dp_rel_;
310  OpmLog::info(os.str());
311  }
312  return resmax < tol_p;
313  }
314 
315  };
316 
317 
318 
319 
320 
322  template <class Grid, class WellModel>
323  struct ModelTraits< BlackoilPressureModel<Grid, WellModel> >
324  {
325  typedef BlackoilState ReservoirState;
329  };
330 
331 } // namespace Opm
332 
333 
334 
335 #endif // OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
BlackoilPressureModel(const typename Base::ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const StandardWells &std_wells, 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: BlackoilPressureModel.hpp:67
A model implementation for the pressure equation in three-phase black oil.
Definition: BlackoilPressureModel.hpp:41
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once per timestep.
Definition: BlackoilPressureModel.hpp:87
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelBase_impl.hpp:758
Definition: AutoDiff.hpp:297
virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual &residual) const =0
Solve the linear system Ax = b, with A being the combined derivative matrix of the residual and b bei...
int sizeNonLinear() const
The size of the non-linear system.
Definition: LinearisedBlackoilResidual.cpp:22
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelBase_impl.hpp:371
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilModelBase_impl.hpp:1148
ADB well_eq
The well_eq has size equal to the number of wells.
Definition: LinearisedBlackoilResidual.hpp:64
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
Class for handling the standard well model.
Definition: StandardWells.hpp:51
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
std::vector< ADB > material_balance_eq
The material_balance_eq vector has one element for each active phase, each of which has size equal to...
Definition: LinearisedBlackoilResidual.hpp:54
Traits to encapsulate the types used by classes using or extending this model.
Definition: BlackoilModelBase.hpp:60
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
ADB well_flux_eq
The well_flux_eq has size equal to the number of wells times the number of phases.
Definition: LinearisedBlackoilResidual.hpp:59
Residual structure of the fully implicit solver.
Definition: LinearisedBlackoilResidual.hpp:47
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
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Struct for containing iteration variables.
Definition: DefaultBlackoilSolutionState.hpp:29
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
BlackoilPressureModel< Grid, WellModel > & asImpl()
Access the most-derived class used for static polymorphism (CRTP).
Definition: BlackoilModelBase.hpp:370
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:383
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilPressureModel.hpp:101
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
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Construct a set of primary variables, each initialized to a given vector.
Definition: AutoDiffBlock.hpp:203
int numMaterials() const
The number of active materials in the model.
Definition: BlackoilModelBase_impl.hpp:394
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415