diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 73dd738ff..4b96ead72 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -45,15 +45,10 @@ namespace Opm const double* porosity, std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - std::transform(porosity, porosity + num_cells, - grid.cell_volumes, - porevol.begin(), - std::multiplies()); + computePorevolume(grid.number_of_cells, grid.cell_volumes, + porosity, porevol); } - /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity array of grid.number_of_cells porosity values @@ -66,13 +61,11 @@ namespace Opm const std::vector& pressure, std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - for (int i = 0; i < num_cells; ++i) { - porevol[i] = porosity[i]*grid.cell_volumes[i]*rock_comp.poroMult(pressure[i]); - } + computePorevolume(grid.number_of_cells, grid.cell_volumes, porosity, rock_comp, pressure, + porevol); } + /// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity_standard array of grid.number_of_cells porosity values (at standard conditions) @@ -482,44 +475,8 @@ namespace Opm const double* densities, const double gravity, const bool per_grid_cell, std::vector& wdp) { - const int nw = wells.number_of_wells; - const size_t np = per_grid_cell ? - saturations.size()/grid.number_of_cells - : saturations.size()/wells.well_connpos[nw]; - // Simple for now: - for (int i = 0; i < nw; i++) { - double depth_ref = wells.depth_ref[i]; - for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) { - int cell = wells.well_cells[j]; - - // Is this correct wrt. depth_ref? - double cell_depth = grid.cell_centroids[3 * cell + 2]; - - double saturation_sum = 0.0; - for (size_t p = 0; p < np; p++) { - if (!per_grid_cell) { - saturation_sum += saturations[j * np + p]; - } else { - saturation_sum += saturations[np * cell + p]; - } - } - if (saturation_sum == 0) { - saturation_sum = 1.0; - } - double density = 0.0; - for (size_t p = 0; p < np; p++) { - if (!per_grid_cell) { - density += saturations[j * np + p] * densities[p] / saturation_sum; - } else { - // Is this a smart way of doing it? - density += saturations[np * cell + p] * densities[p] / saturation_sum; - } - } - - // Is the sign correct? - wdp.push_back(density * (cell_depth - depth_ref) * gravity); - } - } + computeWDP(wells, grid.number_of_cells, grid.cell_centroids, saturations, densities, + gravity, per_grid_cell, wdp); } diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index 1a4a69c79..1362e84d7 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -41,6 +41,16 @@ namespace Opm const double* porosity, std::vector& porevol); + /// @brief Computes pore volume of all cells in a grid. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] begin_cell_volume Iterator to the volume of the first cell. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + std::vector& porevol); /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid @@ -54,6 +64,21 @@ namespace Opm const std::vector& pressure, std::vector& porevol); + /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] Pointer to/ Iterator at the first cell volume. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[in] rock_comp rock compressibility properties + /// @param[in] pressure pressure by cell + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + const RockCompressibility& rock_comp, + const std::vector& pressure, + std::vector& porevol); + /// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @param[in] grid a grid /// @param[in] porosity_standard array of grid.number_of_cells porosity values (at reference presure) @@ -238,6 +263,24 @@ namespace Opm const double* densities, const double gravity, const bool per_grid_cell, std::vector& wdp); + /// Computes the WDP for each well. + /// \param[in] wells Wells that need their wdp calculated. + /// \param[in] number_of_cells The number of cells in the grid. + /// \param[in] begin_cell_centroids Pointer/Iterator to the first cell centroid. + /// \param[in] saturations A vector of weights for each cell for each phase + /// in the grid (or well, see per_grid_cell parameter). So for cell i, + /// saturations[i*densities.size() + p] should give the weight + /// of phase p in cell i. + /// \param[in] densities Density for each phase. + /// \param[out] wdp Will contain, for each well, the wdp of the well. + /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per + /// well cell. + template + void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, + const std::vector& saturations, + const double* densities, const double gravity, const bool per_grid_cell, + std::vector& wdp); + /// Computes (sums) the flow rate for each well. /// \param[in] wells The wells for which the flow rate should be computed. /// \param[in] flow_rates_per_cell Flow rates per well cells. Should ordered the same way as diff --git a/opm/core/utility/miscUtilities_impl.hpp b/opm/core/utility/miscUtilities_impl.hpp index 9f1e693fe..25bc56111 100644 --- a/opm/core/utility/miscUtilities_impl.hpp +++ b/opm/core/utility/miscUtilities_impl.hpp @@ -1,4 +1,7 @@ #include +#include +#include + namespace Opm { /// @brief Estimates a scalar cell velocity from face fluxes. @@ -38,4 +41,83 @@ namespace Opm } } } + + template + void computePorevolume(int number_of_cells, + T begin_cell_volume, + const double* porosity, + std::vector& porevol) + { + porevol.resize(number_of_cells); + std::transform(porosity, porosity + number_of_cells, + begin_cell_volume, + porevol.begin(), + std::multiplies()); + } + + /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. + /// @param[in] number_of_cells The number of cells of the grid. + /// @param[in] porosity array of grid.number_of_cells porosity values + /// @param[in] rock_comp rock compressibility properties + /// @param[in] pressure pressure by cell + /// @param[out] porevol the pore volume by cell. + template + void computePorevolume(int number_of_cells, + T begin_cell_volumes, + const double* porosity, + const RockCompressibility& rock_comp, + const std::vector& pressure, + std::vector& porevol) + { + porevol.resize(number_of_cells); + for (int i = 0; i < number_of_cells; ++i) { + porevol[i] = porosity[i]*begin_cell_volumes[i]*rock_comp.poroMult(pressure[i]); + } + } + + template + void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, const std::vector& saturations, + const double* densities, const double gravity, const bool per_grid_cell, + std::vector& wdp) + { + const int nw = wells.number_of_wells; + const size_t np = per_grid_cell ? + saturations.size()/number_of_cells + : saturations.size()/wells.well_connpos[nw]; + // Simple for now: + for (int i = 0; i < nw; i++) { + double depth_ref = wells.depth_ref[i]; + for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) { + int cell = wells.well_cells[j]; + + // Is this correct wrt. depth_ref? + double cell_depth = UgGridHelpers + ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, cell, 3), 2); + + double saturation_sum = 0.0; + for (size_t p = 0; p < np; p++) { + if (!per_grid_cell) { + saturation_sum += saturations[j * np + p]; + } else { + saturation_sum += saturations[np * cell + p]; + } + } + if (saturation_sum == 0) { + saturation_sum = 1.0; + } + double density = 0.0; + for (size_t p = 0; p < np; p++) { + if (!per_grid_cell) { + density += saturations[j * np + p] * densities[p] / saturation_sum; + } else { + // Is this a smart way of doing it? + density += saturations[np * cell + p] * densities[p] / saturation_sum; + } + } + + // Is the sign correct? + wdp.push_back(density * (cell_depth - depth_ref) * gravity); + } + } + } }