00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED
00021 #define OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED
00022
00023 #include <opm/autodiff/BlackoilMultiSegmentModel.hpp>
00024
00025 #include <opm/autodiff/AutoDiffBlock.hpp>
00026 #include <opm/autodiff/AutoDiffHelpers.hpp>
00027 #include <opm/autodiff/GridHelpers.hpp>
00028 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00029 #include <opm/autodiff/GeoProps.hpp>
00030 #include <opm/autodiff/WellDensitySegmented.hpp>
00031 #include <opm/autodiff/VFPProperties.hpp>
00032 #include <opm/autodiff/VFPProdProperties.hpp>
00033 #include <opm/autodiff/VFPInjProperties.hpp>
00034
00035 #include <opm/core/grid.h>
00036 #include <opm/core/linalg/LinearSolverInterface.hpp>
00037 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00038 #include <opm/core/props/rock/RockCompressibility.hpp>
00039 #include <opm/common/ErrorMacros.hpp>
00040 #include <opm/common/Exceptions.hpp>
00041 #include <opm/parser/eclipse/Units/Units.hpp>
00042 #include <opm/core/well_controls.h>
00043 #include <opm/core/utility/parameters/ParameterGroup.hpp>
00044 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00045
00046 #include <cassert>
00047 #include <cmath>
00048 #include <iostream>
00049 #include <iomanip>
00050 #include <limits>
00051 #include <vector>
00052
00053
00054
00055 namespace Opm {
00056
00057
00058 template <class Grid>
00059 BlackoilMultiSegmentModel<Grid>::
00060 BlackoilMultiSegmentModel(const typename Base::ModelParameters& param,
00061 const Grid& grid ,
00062 const BlackoilPropsAdFromDeck& fluid,
00063 const DerivedGeology& geo ,
00064 const RockCompressibility* rock_comp_props,
00065 const MultisegmentWells& well_model,
00066 const NewtonIterationBlackoilInterface& linsolver,
00067 std::shared_ptr< const EclipseState > eclState,
00068 const bool has_disgas,
00069 const bool has_vapoil,
00070 const bool terminal_output)
00071 : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
00072 eclState, has_disgas, has_vapoil, terminal_output)
00073 {
00074 }
00075
00076
00077
00078
00079
00080 template <class Grid>
00081 void
00082 BlackoilMultiSegmentModel<Grid>::
00083 prepareStep(const SimulatorTimerInterface& timer,
00084 const ReservoirState& reservoir_state,
00085 const WellState& well_state)
00086 {
00087 const double dt = timer.currentStepLength();
00088 pvdt_ = geo_.poreVolume() / dt;
00089 if (active_[Gas]) {
00090 updatePrimalVariableFromState(reservoir_state);
00091 }
00092
00093 const int nw = wellsMultiSegment().size();
00094
00095 if ( !msWellOps().has_multisegment_wells ) {
00096 wellModel().segVDt() = V::Zero(nw);
00097 return;
00098 }
00099
00100 const int nseg_total = well_state.numSegments();
00101 std::vector<double> segment_volume;
00102 segment_volume.reserve(nseg_total);
00103 for (int w = 0; w < nw; ++w) {
00104 WellMultiSegmentConstPtr well = wellsMultiSegment()[w];
00105 const std::vector<double>& segment_volume_well = well->segmentVolume();
00106 segment_volume.insert(segment_volume.end(), segment_volume_well.begin(), segment_volume_well.end());
00107 }
00108 assert(int(segment_volume.size()) == nseg_total);
00109 wellModel().segVDt() = Eigen::Map<V>(segment_volume.data(), nseg_total) / dt;
00110 }
00111
00112
00113
00114
00115
00116 template <class Grid>
00117 void
00118 BlackoilMultiSegmentModel<Grid>::makeConstantState(SolutionState& state) const
00119 {
00120 Base::makeConstantState(state);
00121 state.segp = ADB::constant(state.segp.value());
00122 state.segqs = ADB::constant(state.segqs.value());
00123 }
00124
00125
00126
00127
00128
00129 template <class Grid>
00130 SimulatorReport
00131 BlackoilMultiSegmentModel<Grid>::
00132 assemble(const ReservoirState& reservoir_state,
00133 WellState& well_state,
00134 const bool initial_assembly)
00135 {
00136 using namespace Opm::AutoDiffGrid;
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 wellModel().updateWellControls(well_state);
00152
00153
00154 if (asImpl().wellModel().wellCollection()->groupControlActive()) {
00155
00156 Base::applyVREPGroupControl(reservoir_state, well_state);
00157
00158 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
00159 }
00160
00161
00162 SolutionState state = asImpl().variableState(reservoir_state, well_state);
00163
00164 if (initial_assembly) {
00165
00166 SolutionState state0 = state;
00167 asImpl().makeConstantState(state0);
00168
00169
00170 asImpl().computeAccum(state0, 0);
00171 wellModel().computeSegmentFluidProperties(state0);
00172 const int np = numPhases();
00173 assert(np == int(wellModel().segmentCompSurfVolumeInitial().size()));
00174 for (int phase = 0; phase < np; ++phase) {
00175 wellModel().segmentCompSurfVolumeInitial()[phase] = wellModel().segmentCompSurfVolumeCurrent()[phase].value();
00176 }
00177
00178 asImpl().computeWellConnectionPressures(state0, well_state);
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 asImpl().assembleMassBalanceEq(state);
00192
00193
00194 if ( ! wellsActive() ) {
00195 SimulatorReport report;
00196 return report;
00197 }
00198
00199 wellModel().computeSegmentFluidProperties(state);
00200
00201 const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
00202 wellModel().computeSegmentPressuresDelta(gravity);
00203
00204 std::vector<ADB> mob_perfcells;
00205 std::vector<ADB> b_perfcells;
00206 SimulatorReport report;
00207 wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
00208 if (param_.solve_welleq_initially_ && initial_assembly) {
00209
00210 report = asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
00211 }
00212
00213
00214
00215 V aliveWells;
00216 std::vector<ADB> cq_s;
00217 wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
00218 wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
00219 wellModel().addWellFluxEq(cq_s, state, residual_);
00220 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
00221 wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
00222 return report;
00223 }
00224
00225
00226
00227
00228
00229 template <class Grid>
00230 SimulatorReport
00231 BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
00232 const std::vector<ADB>& b_perfcells,
00233 const ReservoirState& reservoir_state,
00234 SolutionState& state,
00235 WellState& well_state)
00236 {
00237 SimulatorReport report = Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
00238
00239 if (report.converged) {
00240
00241
00242 const int np = numPhases();
00243 const int nseg_total =well_state.numSegments();
00244 {
00245
00246
00247 ADB::V new_segp = Eigen::Map<ADB::V>(well_state.segPress().data(), nseg_total);
00248
00249
00250 std::vector<ADB::M> old_segp_derivs = state.segp.derivative();
00251 state.segp = ADB::function(std::move(new_segp), std::move(old_segp_derivs));
00252 }
00253 {
00254
00255
00256
00257 const DataBlock segrates = Eigen::Map<const DataBlock>(well_state.segPhaseRates().data(), nseg_total, np).transpose();
00258 ADB::V new_segqs = Eigen::Map<const V>(segrates.data(), nseg_total * np);
00259 std::vector<ADB::M> old_segqs_derivs = state.segqs.derivative();
00260 state.segqs = ADB::function(std::move(new_segqs), std::move(old_segqs_derivs));
00261 }
00262
00263
00264
00265 asImpl().computeWellConnectionPressures(state, well_state);
00266 }
00267
00268 return report;
00269 }
00270
00271
00272
00273
00274 template <class Grid>
00275 void
00276 BlackoilMultiSegmentModel<Grid>::
00277 computeWellConnectionPressures(const SolutionState& state,
00278 const WellState& well_state)
00279 {
00280 const int np = numPhases();
00281 const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
00282 std::vector<ADB> fluid_density(np, ADB::null());
00283
00284 for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
00285 const int canonicalPhaseIdx = canph_[phaseIdx];
00286 fluid_density[phaseIdx] = fluidDensity(canonicalPhaseIdx, sd_.rq[phaseIdx].b, state.rs, state.rv);
00287 }
00288 wellModel().computeWellConnectionPressures(state, well_state, kr_adb, fluid_density);
00289 }
00290
00291 }
00292
00293 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED