00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
00022 #define OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
00023
00024
00025 #include <opm/autodiff/BlackoilModelBase.hpp>
00026 #include <opm/core/simulator/BlackoilState.hpp>
00027 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00028 #include <opm/autodiff/BlackoilModelParameters.hpp>
00029 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
00030
00031 #include <algorithm>
00032
00033 namespace Opm {
00034
00040 template<class Grid, class WellModel>
00041 class BlackoilPressureModel : public BlackoilModelBase<Grid, WellModel, BlackoilPressureModel<Grid, WellModel> >
00042 {
00043 public:
00044
00045 typedef BlackoilModelBase<Grid, WellModel, BlackoilPressureModel<Grid, WellModel> > Base;
00046 friend Base;
00047
00048 typedef typename Base::ReservoirState ReservoirState;
00049 typedef typename Base::WellState WellState;
00050 typedef typename Base::SolutionState SolutionState;
00051 typedef typename Base::V V;
00052
00067 BlackoilPressureModel(const typename Base::ModelParameters& param,
00068 const Grid& grid,
00069 const BlackoilPropsAdFromDeck& fluid,
00070 const DerivedGeology& geo,
00071 const RockCompressibility* rock_comp_props,
00072 const StandardWells& std_wells,
00073 const NewtonIterationBlackoilInterface& linsolver,
00074 std::shared_ptr< const EclipseState> eclState,
00075 const bool has_disgas,
00076 const bool has_vapoil,
00077 const bool terminal_output)
00078 : Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
00079 eclState, has_disgas, has_vapoil, terminal_output),
00080 state0_(3),
00081 max_dp_rel_(std::numeric_limits<double>::infinity()),
00082 scaling_{ ADB::null(), ADB::null(), ADB::null() }
00083 {
00084 }
00085
00087 void prepareStep(const SimulatorTimerInterface& timer,
00088 const ReservoirState& reservoir_state,
00089 const WellState& well_state)
00090 {
00091 asImpl().wellModel().setStoreWellPerforationFluxesFlag(true);
00092 Base::prepareStep(timer, reservoir_state, well_state);
00093 max_dp_rel_ = std::numeric_limits<double>::infinity();
00094 state0_ = asImpl().variableState(reservoir_state, well_state);
00095 asImpl().makeConstantState(state0_);
00096 }
00097
00098
00101 V solveJacobianSystem() const
00102 {
00103
00104
00105
00106 const int n1 = residual_.material_balance_eq[0].size();
00107 const int n2 = residual_.well_flux_eq.size() + residual_.well_eq.size();
00108 const int n_full = residual_.sizeNonLinear();
00109 const auto& mb = residual_.material_balance_eq;
00110 const auto& fe = residual_.well_flux_eq;
00111 const auto& we = residual_.well_eq;
00112 LinearisedBlackoilResidual pressure_res = {
00113 {
00114
00115 ADB::function(mb[0].value(), { mb[0].derivative()[0], mb[0].derivative()[3], mb[0].derivative()[4] })
00116 },
00117 ADB::function(fe.value(), { fe.derivative()[0], fe.derivative()[3], fe.derivative()[4] }),
00118 ADB::function(we.value(), { we.derivative()[0], we.derivative()[3], we.derivative()[4] }),
00119 residual_.matbalscale,
00120 residual_.singlePrecision
00121 };
00122 assert(pressure_res.sizeNonLinear() == n1 + n2);
00123 V dx_pressure = linsolver_.computeNewtonIncrement(pressure_res);
00124 assert(dx_pressure.size() == n1 + n2);
00125 V dx_full = V::Zero(n_full);
00126 dx_full.topRows(n1) = dx_pressure.topRows(n1);
00127 dx_full.bottomRows(n2) = dx_pressure.bottomRows(n2);
00128 return dx_full;
00129 }
00130
00131 using Base::numPhases;
00132 using Base::numMaterials;
00133 using Base::wellModel;
00134
00135 protected:
00136 using Base::asImpl;
00137 using Base::linsolver_;
00138 using Base::residual_;
00139 using Base::sd_;
00140 using Base::grid_;
00141 using Base::ops_;
00142 using Base::has_vapoil_;
00143 using Base::has_disgas_;
00144
00145 SolutionState state0_;
00146 double max_dp_rel_ = std::numeric_limits<double>::infinity();
00147 ADB scaling_[3] = { ADB::null(), ADB::null(), ADB::null() };
00148
00149
00150
00151
00152 SimulatorReport
00153 assemble(const ReservoirState& reservoir_state,
00154 WellState& well_state,
00155 const bool initial_assembly)
00156 {
00157 SimulatorReport report;
00158
00159 report += Base::assemble(reservoir_state, well_state, initial_assembly);
00160 if (initial_assembly) {
00161 }
00162
00163
00164 ADB pressure_residual = ADB::constant(V::Zero(residual_.material_balance_eq[0].size()));
00165 for (int phase = 0; phase < numPhases(); ++phase) {
00166 pressure_residual += residual_.material_balance_eq[phase] * scaling_[phase];
00167 }
00168 residual_.material_balance_eq[0] = pressure_residual;
00169
00170
00171 const int n = sd_.rq[0].mflux.size();
00172 V flux = V::Zero(n);
00173 for (int phase = 0; phase < numPhases(); ++phase) {
00174 UpwindSelector<double> upwind(grid_, ops_, sd_.rq[phase].dh.value());
00175 flux += sd_.rq[phase].mflux.value() / upwind.select(sd_.rq[phase].b.value());
00176 }
00177
00178
00179
00180
00181
00182 ReservoirState& s = const_cast<ReservoirState&>(reservoir_state);
00183 s.faceflux().resize(n);
00184 std::copy_n(flux.data(), n, s.faceflux().begin());
00185 if (asImpl().localWellsActive()) {
00186 const V& wflux = asImpl().wellModel().getStoredWellPerforationFluxes();
00187 assert(int(well_state.perfRates().size()) == wflux.size());
00188 std::copy_n(wflux.data(), wflux.size(), well_state.perfRates().begin());
00189 }
00190
00191 return report;
00192 }
00193
00194
00195
00196
00197
00198 SolutionState
00199 variableState(const ReservoirState& x,
00200 const WellState& xw) const
00201 {
00202
00203 std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
00204 std::vector<ADB> vars = ADB::variables(vars0);
00205 const std::vector<int> indices = asImpl().variableStateIndices();
00206 vars[indices[Sw]] = ADB::constant(vars[indices[Sw]].value());
00207 vars[indices[Xvar]] = ADB::constant(vars[indices[Xvar]].value());
00208
00209 return asImpl().variableStateExtractVars(x, indices, vars);
00210 }
00211
00212
00213
00214
00215
00216 void computeAccum(const SolutionState& state,
00217 const int aix )
00218 {
00219 if (aix == 0) {
00220 Base::computeAccum(state0_, aix);
00221 } else {
00222 Base::computeAccum(state, aix);
00223 }
00224 }
00225
00226
00227
00228
00229
00230 void assembleMassBalanceEq(const SolutionState& state)
00231 {
00232 Base::assembleMassBalanceEq(state);
00233
00234
00235
00236
00237 assert(numPhases() == 3);
00238 assert(numMaterials() == 3);
00239 V one = V::Constant(state.pressure.size(), 1.0);
00240 scaling_[Water] = one / sd_.rq[Water].b;
00241 scaling_[Oil] = one / sd_.rq[Oil].b - state.rs / sd_.rq[Gas].b;
00242 scaling_[Gas] = one / sd_.rq[Gas].b - state.rv / sd_.rq[Oil].b;
00243 if (has_disgas_ && has_vapoil_) {
00244 ADB r_factor = one / (one - state.rs * state.rv);
00245 scaling_[Oil] = r_factor * scaling_[Oil];
00246 scaling_[Gas] = r_factor * scaling_[Gas];
00247 }
00248
00249
00250 }
00251
00252
00253
00254
00255
00256 void updateState(const V& dx,
00257 ReservoirState& reservoir_state,
00258 WellState& well_state)
00259 {
00260
00261
00262 std::vector<double> rs_old = reservoir_state.gasoilratio();
00263 std::vector<double> rv_old = reservoir_state.rv();
00264 auto hs_old = reservoir_state.hydroCarbonState();
00265 auto phasecond_old = Base::phaseCondition_;
00266 auto isRs_old = Base::isRs_;
00267 auto isRv_old = Base::isRv_;
00268 auto isSg_old = Base::isSg_;
00269
00270
00271 const auto minmax_iters = std::minmax_element(reservoir_state.pressure().begin(),
00272 reservoir_state.pressure().end());
00273 const double range = *minmax_iters.second - *minmax_iters.first;
00274
00275
00276 Base::updateState(dx, reservoir_state, well_state);
00277
00278
00279 max_dp_rel_ = dx.head(reservoir_state.pressure().size()).abs().maxCoeff() / range;
00280
00281
00282 reservoir_state.gasoilratio() = rs_old;
00283 reservoir_state.rv() = rv_old;
00284 reservoir_state.hydroCarbonState() = hs_old;
00285 Base::phaseCondition_ = phasecond_old;
00286 Base::isRs_ = isRs_old;
00287 Base::isRv_ = isRv_old;
00288 Base::isSg_ = isSg_old;
00289 }
00290
00291
00292
00293
00294
00295 bool getConvergence(const SimulatorTimerInterface& , const int iteration)
00296 {
00297 const double tol_p = 1e-10;
00298 const double resmax = residual_.material_balance_eq[0].value().abs().maxCoeff();
00299 if (Base::terminalOutputEnabled()) {
00300
00301 if (iteration == 0) {
00302 OpmLog::info("\nIter Res(p) Delta(p)\n");
00303 }
00304 std::ostringstream os;
00305 os.precision(3);
00306 os.setf(std::ios::scientific);
00307 os << std::setw(4) << iteration;
00308 os << std::setw(11) << resmax;
00309 os << std::setw(11) << max_dp_rel_;
00310 OpmLog::info(os.str());
00311 }
00312 return resmax < tol_p;
00313 }
00314
00315 };
00316
00317
00318
00319
00320
00322 template <class Grid, class WellModel>
00323 struct ModelTraits< BlackoilPressureModel<Grid, WellModel> >
00324 {
00325 typedef BlackoilState ReservoirState;
00326 typedef WellStateFullyImplicitBlackoil WellState;
00327 typedef BlackoilModelParameters ModelParameters;
00328 typedef DefaultBlackoilSolutionState SolutionState;
00329 };
00330
00331 }
00332
00333
00334
00335 #endif // OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED