Added support for computePoreVolume for grids apart from UnstructuredGrid.

This commit is contained in:
Markus Blatt 2014-02-25 15:12:16 +01:00
parent e30031dc77
commit ab0b2a9c6b
3 changed files with 132 additions and 50 deletions

View File

@ -45,15 +45,10 @@ namespace Opm
const double* porosity, const double* porosity,
std::vector<double>& porevol) std::vector<double>& porevol)
{ {
int num_cells = grid.number_of_cells; computePorevolume(grid.number_of_cells, grid.cell_volumes,
porevol.resize(num_cells); porosity, porevol);
std::transform(porosity, porosity + num_cells,
grid.cell_volumes,
porevol.begin(),
std::multiplies<double>());
} }
/// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
/// @param[in] grid a grid /// @param[in] grid a grid
/// @param[in] porosity array of grid.number_of_cells porosity values /// @param[in] porosity array of grid.number_of_cells porosity values
@ -66,13 +61,11 @@ namespace Opm
const std::vector<double>& pressure, const std::vector<double>& pressure,
std::vector<double>& porevol) std::vector<double>& porevol)
{ {
int num_cells = grid.number_of_cells; computePorevolume(grid.number_of_cells, grid.cell_volumes, porosity, rock_comp, pressure,
porevol.resize(num_cells); porevol);
for (int i = 0; i < num_cells; ++i) {
porevol[i] = porosity[i]*grid.cell_volumes[i]*rock_comp.poroMult(pressure[i]);
}
} }
/// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
/// @param[in] grid a grid /// @param[in] grid a grid
/// @param[in] porosity_standard array of grid.number_of_cells porosity values (at standard conditions) /// @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, const double* densities, const double gravity, const bool per_grid_cell,
std::vector<double>& wdp) std::vector<double>& wdp)
{ {
const int nw = wells.number_of_wells; computeWDP(wells, grid.number_of_cells, grid.cell_centroids, saturations, densities,
const size_t np = per_grid_cell ? gravity, per_grid_cell, wdp);
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);
}
}
} }

View File

@ -41,6 +41,16 @@ namespace Opm
const double* porosity, const double* porosity,
std::vector<double>& porevol); std::vector<double>& 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<class T>
void computePorevolume(int number_of_cells,
T begin_cell_volume,
const double* porosity,
std::vector<double>& porevol);
/// @brief Computes pore volume of all cells in a grid, with rock compressibility effects. /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
/// @param[in] grid a grid /// @param[in] grid a grid
@ -54,6 +64,21 @@ namespace Opm
const std::vector<double>& pressure, const std::vector<double>& pressure,
std::vector<double>& porevol); std::vector<double>& 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<class T>
void computePorevolume(int number_of_cells,
T begin_cell_volume,
const double* porosity,
const RockCompressibility& rock_comp,
const std::vector<double>& pressure,
std::vector<double>& porevol);
/// @brief Computes porosity of all cells in a grid, with rock compressibility effects. /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
/// @param[in] grid a grid /// @param[in] grid a grid
/// @param[in] porosity_standard array of grid.number_of_cells porosity values (at reference presure) /// @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, const double* densities, const double gravity, const bool per_grid_cell,
std::vector<double>& wdp); std::vector<double>& 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<class T>
void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids,
const std::vector<double>& saturations,
const double* densities, const double gravity, const bool per_grid_cell,
std::vector<double>& wdp);
/// Computes (sums) the flow rate for each well. /// Computes (sums) the flow rate for each well.
/// \param[in] wells The wells for which the flow rate should be computed. /// \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 /// \param[in] flow_rates_per_cell Flow rates per well cells. Should ordered the same way as

View File

@ -1,4 +1,7 @@
#include <opm/core/grid/GridHelpers.hpp> #include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/wells.h>
#include <opm/core/props/rock/RockCompressibility.hpp>
namespace Opm namespace Opm
{ {
/// @brief Estimates a scalar cell velocity from face fluxes. /// @brief Estimates a scalar cell velocity from face fluxes.
@ -38,4 +41,83 @@ namespace Opm
} }
} }
} }
template<class T>
void computePorevolume(int number_of_cells,
T begin_cell_volume,
const double* porosity,
std::vector<double>& porevol)
{
porevol.resize(number_of_cells);
std::transform(porosity, porosity + number_of_cells,
begin_cell_volume,
porevol.begin(),
std::multiplies<double>());
}
/// @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<class T>
void computePorevolume(int number_of_cells,
T begin_cell_volumes,
const double* porosity,
const RockCompressibility& rock_comp,
const std::vector<double>& pressure,
std::vector<double>& 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<class T>
void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, const std::vector<double>& saturations,
const double* densities, const double gravity, const bool per_grid_cell,
std::vector<double>& 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);
}
}
}
} }