24 template <
class Gr
idT>
25 SimulatorFullyImplicitBlackoilPolymer<GridT>::
26 SimulatorFullyImplicitBlackoilPolymer(
const ParameterGroup& param,
29 BlackoilPropsAdFromDeck& props,
30 const PolymerPropsAd& polymer_props,
31 const RockCompressibility* rock_comp_props,
32 NewtonIterationBlackoilInterface& linsolver,
33 const double* gravity,
34 const bool has_disgas,
35 const bool has_vapoil,
36 const bool has_polymer,
37 const bool has_plyshlog,
38 const bool has_shrate,
39 std::shared_ptr<EclipseState> eclipse_state,
40 BlackoilOutputWriter& output_writer,
41 std::shared_ptr< Deck > deck,
42 const std::vector<double>& threshold_pressures_by_face)
54 threshold_pressures_by_face,
56 std::unordered_set<
std::string>())
57 , polymer_props_(polymer_props)
58 , has_polymer_(has_polymer)
59 , has_plyshlog_(has_plyshlog)
60 , has_shrate_(has_shrate)
65 template <
class Gr
idT>
66 auto SimulatorFullyImplicitBlackoilPolymer<GridT>::
67 createSolver(
const WellModel& well_model)
68 -> std::unique_ptr<Solver>
70 typedef typename Traits::Model Model;
73 auto model = std::unique_ptr<Model>(
new Model(BaseType::model_param_,
77 BaseType::rock_comp_props_,
81 BaseType::eclipse_state_,
82 BaseType::has_disgas_,
83 BaseType::has_vapoil_,
90 BaseType::terminal_output_));
92 if (!BaseType::threshold_pressures_by_face_.empty()) {
93 model->setThresholdPressures(BaseType::threshold_pressures_by_face_);
96 return std::unique_ptr<Solver>(
new Solver(BaseType::solver_param_, std::move(model)));
102 template <
class Gr
idT>
103 void SimulatorFullyImplicitBlackoilPolymer<GridT>::
104 handleAdditionalWellInflow(SimulatorTimer& timer,
105 WellsManager& wells_manager,
106 typename BaseType::WellState& well_state,
110 std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
111 if (deck_->hasKeyword(
"WPOLYMER")) {
112 if (wells_manager.c_wells() == 0) {
113 OPM_THROW(std::runtime_error,
"Cannot control polymer injection via WPOLYMER without wells.");
115 polymer_inflow_ptr.reset(
new PolymerInflowFromDeck(*BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
117 OPM_MESSAGE(
"Warning: simulating with no WPOLYMER in deck (no polymer will be injected).");
118 polymer_inflow_ptr.reset(
new PolymerInflowBasic(0.0*Opm::unit::day,
122 std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
123 polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
124 timer.simulationTimeElapsed() + timer.currentStepLength(),
126 well_state.polymerInflow() = polymer_inflow_c;
129 computeRepRadiusPerfLength(*BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_, wells_bore_diameter_);
134 template <
class Gr
idT>
135 void SimulatorFullyImplicitBlackoilPolymer<GridT>::
136 setupCompressedToCartesian(
const int* global_cell,
int number_of_cells,
137 std::map<int,int>& cartesian_to_compressed )
140 for (
int i = 0; i < number_of_cells; ++i) {
141 cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
145 for (
int i = 0; i < number_of_cells; ++i) {
146 cartesian_to_compressed.insert(std::make_pair(i, i));
153 template <
class Gr
idT>
154 void SimulatorFullyImplicitBlackoilPolymer<GridT>::
155 computeRepRadiusPerfLength(
const Opm::EclipseState& eclipseState,
156 const size_t timeStep,
158 std::vector<double>& wells_rep_radius,
159 std::vector<double>& wells_perf_length,
160 std::vector<double>& wells_bore_diameter)
165 int number_of_cells = Opm::UgGridHelpers::numCells(grid);
166 const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
167 const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
168 auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
169 auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
171 if (eclipseState.getSchedule().numWells() == 0) {
172 OPM_MESSAGE(
"No wells specified in Schedule section, " 173 "initializing no wells");
177 const size_t n_perf = wells_rep_radius.size();
179 wells_rep_radius.clear();
180 wells_perf_length.clear();
181 wells_bore_diameter.clear();
183 wells_rep_radius.reserve(n_perf);
184 wells_perf_length.reserve(n_perf);
185 wells_bore_diameter.reserve(n_perf);
187 std::map<int,int> cartesian_to_compressed;
189 setupCompressedToCartesian(global_cell, number_of_cells,
190 cartesian_to_compressed);
192 const auto& schedule = eclipseState.getSchedule();
193 auto wells = schedule.getWells(timeStep);
197 for (
auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
198 const auto* well = (*wellIter);
200 if (well->getStatus(timeStep) == WellCommon::SHUT) {
204 const auto& completionSet = well->getCompletions(timeStep);
205 for (
size_t c=0; c<completionSet.size(); c++) {
206 const auto& completion = completionSet.get(c);
207 if (completion.getState() == WellCompletion::OPEN) {
208 int i = completion.getI();
209 int j = completion.getJ();
210 int k = completion.getK();
212 const int* cpgdim = cart_dims;
213 int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
214 std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
215 if (cgit == cartesian_to_compressed.end()) {
216 OPM_THROW(std::runtime_error,
"Cell with i,j,k indices " << i <<
' ' << j <<
' ' 217 << k <<
" not found in grid (well = " << well->name() <<
')');
219 int cell = cgit->second;
222 double radius = 0.5*completion.getDiameter();
224 radius = 0.5*unit::feet;
225 OPM_MESSAGE(
"**** Warning: Well bore internal radius set to " << radius);
228 const std::array<double, 3> cubical =
229 WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
231 WellCompletion::DirectionEnum direction = completion.getDirection();
237 case Opm::WellCompletion::DirectionEnum::X:
238 re = std::sqrt(cubical[1] * cubical[2] / M_PI);
239 perf_length = cubical[0];
241 case Opm::WellCompletion::DirectionEnum::Y:
242 re = std::sqrt(cubical[0] * cubical[2] / M_PI);
243 perf_length = cubical[1];
245 case Opm::WellCompletion::DirectionEnum::Z:
246 re = std::sqrt(cubical[0] * cubical[1] / M_PI);
247 perf_length = cubical[2];
250 OPM_THROW(std::runtime_error,
" Dirtecion of well is not supported ");
253 double repR = std::sqrt(re * radius);
254 wells_rep_radius.push_back(repR);
255 wells_perf_length.push_back(perf_length);
256 wells_bore_diameter.push_back(2. * radius);
259 if (completion.getState() != WellCompletion::SHUT) {
260 OPM_THROW(std::runtime_error,
"Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) <<
" not handled");
Definition: AutoDiff.hpp:297
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22