From 8205ad3303275377c796393f820447b4b7a95e6e Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 2 Jun 2015 11:09:56 +0200 Subject: [PATCH] calculating the representative radius and perf length for all well perforations, to be used in shear-thinning calculation. The calculation is approximated. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 10 ++ .../BlackoilPolymerModel_impl.hpp | 12 +- .../SimulatorFullyImplicitBlackoilPolymer.hpp | 22 ++- ...latorFullyImplicitBlackoilPolymer_impl.hpp | 134 ++++++++++++++++++ 4 files changed, 170 insertions(+), 8 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 1887e6ac5..937e17a18 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -80,6 +80,8 @@ namespace Opm { const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, + std::vector& wells_rep_radius, + std::vector& wells_perf_length, const bool terminal_output); /// Called once before each time step. @@ -128,6 +130,11 @@ namespace Opm { const int poly_pos_; V cmax_; + // representative radius and perforation length of well perforations + // to be used in shear-thinning computation. + std::vector wells_rep_radius_; + std::vector wells_perf_length_; + // Need to declare Base members we want to use here. using Base::grid_; using Base::fluid_; @@ -252,6 +259,9 @@ namespace Opm { std::vector& maxNormWell, int nc, int nw) const; + + + }; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 35dcc3431..8fd123985 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -86,13 +86,17 @@ namespace Opm { const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, + std::vector& wells_rep_radius, + std::vector& wells_perf_length, const bool terminal_output) : Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver, has_disgas, has_vapoil, terminal_output), polymer_props_ad_(polymer_props_ad), has_polymer_(has_polymer), has_plyshlog_(has_plyshlog), - poly_pos_(detail::polymerPos(fluid.phaseUsage())) + poly_pos_(detail::polymerPos(fluid.phaseUsage())), + wells_rep_radius_(wells_rep_radius), + wells_perf_length_(wells_perf_length) { if (has_polymer_) { if (!active_[Water]) { @@ -576,10 +580,6 @@ namespace Opm { return converged; } - - - - template ADB BlackoilPolymerModel::computeMc(const SolutionState& state) const @@ -587,8 +587,6 @@ namespace Opm { return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration); } - - } // namespace Opm #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 7bc71b15f..600dd9768 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -124,17 +124,37 @@ namespace Opm std::unique_ptr createSolver(const Wells* wells); + void handleAdditionalWellInflow(SimulatorTimer& timer, WellsManager& wells_manager, typename BaseType::WellState& well_state, const Wells* wells); + private: const PolymerPropsAd& polymer_props_; bool has_polymer_; - /// flag for PLYSHLOG keywords + // flag for PLYSHLOG keyword bool has_plyshlog_; DeckConstPtr deck_; + + std::vector wells_rep_radius_; + std::vector wells_perf_length_; + + // generate the mapping from Cartesian grid cells to global compressed cells, + // copied from opm-core, to be used in function computeRepRadiusPerfLength() + static void + setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed); + + // calculate the representative radius and length for for well peforations + // it will be used in the shear-thinning calcluation only. + void + computeRepRadiusPerfLength(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, + const GridT& grid, + std::vector& wells_rep_radius, + std::vector& wells_perf_length); + }; } // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 94be35d02..0ef805e4b 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -68,6 +68,7 @@ namespace Opm ModelParams modelParams( BaseType::param_ ); typedef NewtonSolver Solver; + auto model = std::unique_ptr(new Model(modelParams, BaseType::grid_, BaseType::props_, @@ -80,6 +81,8 @@ namespace Opm BaseType::has_vapoil_, has_polymer_, has_plyshlog_, + wells_rep_radius_, + wells_perf_length_, BaseType::terminal_output_)); if (!BaseType::threshold_pressures_by_face_.empty()) { @@ -91,6 +94,9 @@ namespace Opm return std::unique_ptr(new Solver(solverParams, std::move(model))); } + + + template void SimulatorFullyImplicitBlackoilPolymer:: handleAdditionalWellInflow(SimulatorTimer& timer, @@ -115,6 +121,134 @@ namespace Opm timer.simulationTimeElapsed() + timer.currentStepLength(), polymer_inflow_c); well_state.polymerInflow() = polymer_inflow_c; + + computeRepRadiusPerfLength(BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_); + } + + + template + void SimulatorFullyImplicitBlackoilPolymer:: + setupCompressedToCartesian(const int* global_cell, int number_of_cells, + std::map& cartesian_to_compressed ) + { + if (global_cell) { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + } + + + template + void SimulatorFullyImplicitBlackoilPolymer:: + computeRepRadiusPerfLength(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, + const GridT& grid, + std::vector& wells_rep_radius, + std::vector& wells_perf_length) + { + + int number_of_cells = Opm::UgGridHelpers::numCells(grid); + const int* global_cell = Opm::UgGridHelpers::globalCell(grid); + const int* cart_dims = Opm::UgGridHelpers::cartDims(grid); + auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid); + auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid); + + if (eclipseState->getSchedule()->numWells() == 0) { + OPM_MESSAGE("No wells specified in Schedule section, " + "initializing no wells"); + return; + } + + wells_rep_radius.clear(); + wells_perf_length.clear(); + + std::map cartesian_to_compressed; + + setupCompressedToCartesian(global_cell, number_of_cells, + cartesian_to_compressed); + + ScheduleConstPtr schedule = eclipseState->getSchedule(); + std::vector wells = schedule->getWells(timeStep); + + int well_index = 0; + + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + + if (well->getStatus(timeStep) == WellCommon::SHUT) { + continue; + } + { // COMPDAT handling + CompletionSetConstPtr completionSet = well->getCompletions(timeStep); + for (size_t c=0; csize(); c++) { + CompletionConstPtr completion = completionSet->get(c); + if (completion->getState() == WellCompletion::OPEN) { + int i = completion->getI(); + int j = completion->getJ(); + int k = completion->getK(); + + const int* cpgdim = cart_dims; + int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); + std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + } + int cell = cgit->second; + + { + double radius = 0.5*completion->getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + + const std::array cubical = + WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell); + + WellCompletion::DirectionEnum direction = completion->getDirection(); + + double re; // area equivalent radius of the grid block + double perf_length; // the length of the well perforation + + switch (direction) { + case Opm::WellCompletion::DirectionEnum::X: + re = std::sqrt(cubical[1] * cubical[2] / M_PI); + perf_length = cubical[0]; + break; + case Opm::WellCompletion::DirectionEnum::Y: + re = std::sqrt(cubical[0] * cubical[2] / M_PI); + perf_length = cubical[1]; + break; + case Opm::WellCompletion::DirectionEnum::Z: + re = std::sqrt(cubical[0] * cubical[1] / M_PI); + perf_length = cubical[2]; + break; + default: + OPM_THROW(std::runtime_error, " Dirtecion of well is not supported "); + } + + double repR = std::sqrt(re * radius); + wells_rep_radius.push_back(repR); + wells_perf_length.push_back(perf_length); + } + } else { + if (completion->getState() != WellCompletion::SHUT) { + OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion->getState() ) << " not handled"); + } + } + + } + } + well_index++; + } } } // namespace Opm