21 #ifndef OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
22 #define OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
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>
40 template<
class Gr
id,
class WellModel>
48 typedef typename Base::ReservoirState ReservoirState;
49 typedef typename Base::WellState WellState;
50 typedef typename Base::SolutionState SolutionState;
51 typedef typename Base::V V;
71 const RockCompressibility* rock_comp_props,
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),
81 max_dp_rel_(std::numeric_limits<double>::infinity()),
88 const ReservoirState& reservoir_state,
89 const WellState& well_state)
91 asImpl().wellModel().setStoreWellPerforationFluxesFlag(
true);
93 max_dp_rel_ = std::numeric_limits<double>::infinity();
94 state0_ =
asImpl().variableState(reservoir_state, well_state);
95 asImpl().makeConstantState(state0_);
111 const auto& we = residual_.
well_eq;
115 ADB::function(mb[0].value(), { mb[0].derivative()[0], mb[0].derivative()[3], mb[0].derivative()[4] })
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
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);
137 using Base::linsolver_;
138 using Base::residual_;
142 using Base::has_vapoil_;
143 using Base::has_disgas_;
145 SolutionState state0_;
146 double max_dp_rel_ = std::numeric_limits<double>::infinity();
153 assemble(
const ReservoirState& reservoir_state,
154 WellState& well_state,
155 const bool initial_assembly)
157 SimulatorReport report;
159 report +=
Base::assemble(reservoir_state, well_state, initial_assembly);
160 if (initial_assembly) {
165 for (
int phase = 0; phase <
numPhases(); ++phase) {
171 const int n = sd_.rq[0].mflux.size();
173 for (
int phase = 0; phase <
numPhases(); ++phase) {
175 flux += sd_.rq[phase].mflux.value() / upwind.select(sd_.rq[phase].b.value());
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());
199 variableState(
const ReservoirState& x,
200 const WellState& xw)
const
203 std::vector<V> vars0 =
asImpl().variableStateInitials(x, xw);
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());
209 return asImpl().variableStateExtractVars(x, indices, vars);
216 void computeAccum(
const SolutionState& state,
220 Base::computeAccum(state0_, aix);
222 Base::computeAccum(state, aix);
230 void assembleMassBalanceEq(
const SolutionState& state)
232 Base::assembleMassBalanceEq(state);
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];
256 void updateState(
const V& dx,
257 ReservoirState& reservoir_state,
258 WellState& well_state)
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_;
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;
279 max_dp_rel_ = dx.head(reservoir_state.pressure().size()).abs().maxCoeff() / range;
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;
297 const double tol_p = 1e-10;
301 if (iteration == 0) {
302 OpmLog::info(
"\nIter Res(p) Delta(p)\n");
304 std::ostringstream os;
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());
312 return resmax < tol_p;
322 template <
class Gr
id,
class WellModel>
325 typedef BlackoilState ReservoirState;
335 #endif // OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
BlackoilPressureModel(const typename Base::ModelParameters ¶m, 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.
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
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 size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
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.
int sizeNonLinear() const
The size of the non-linear system.
Definition: LinearisedBlackoilResidual.cpp:22
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilPressureModel.hpp:101
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
int numPhases() const
The number of active fluid phases in the model.
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
int numMaterials() const
The number of active materials in the model.
Struct for containing iteration variables.
Definition: DefaultBlackoilSolutionState.hpp:29
Upwind selection in absence of counter-current flow (i.e., without effects of gravity and/or capillar...
Definition: AutoDiffHelpers.hpp:181
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
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
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