OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #89 from atgeirr/gravity-in-wells

Browse Source 
Gravity in wells
This commit is contained in:
Bård Skaflestad 2012-11-07 03:53:04 -08:00
commit d399b2d7cc
7 changed files with 588 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

212
opm/core/newwells.h
View File

@ -33,97 +33,177 @@
extern "C" {
#endif
/** Well type indicates desired/expected well behaviour. */
enum WellType { INJECTOR, PRODUCER };
/** Type of well control equation or inequality constraint.
* BHP -> Well constrained by bottom-hole pressure target.
* RESERVOIR_RATE -> Well constrained by reservoir volume flow rates.
* SURFACE_RATE -> Well constrained by surface volume flow rates.
/**
* Well type indicates desired/expected well behaviour.
*/
enum WellControlType { BHP, RESERVOIR_RATE, SURFACE_RATE };
enum WellType {
INJECTOR, /**< Well is an injector */
PRODUCER /**< Well is a producer */
};
/** Controls for a single well.
* Each control specifies a well rate or bottom-hole pressure. Only
* one control can be active at a time, indicated by current. The
* meaning of each control's target value depends on the control type:
* BHP -> target pressure in Pascal.
* RESERVOIR_RATE -> target reservoir volume rate in cubic(meter)/second
* SURFACE_RATE -> target surface volume rate in cubic(meter)/second
* The sign convention for RATE targets is as follows:
* (+) Fluid flowing into reservoir, i.e. injecting.
* (-) Fluid flowing out of reservoir, i.e. producing.
* For *_RATE controls, the distribution of phases used for the control
* is also needed. For example, a total rate control should have 1.0
* for each phase, whereas a control on oil rate should have 1.0 for
* the oil phase and 0.0 for the rest. For BHP controls, this is unused.
* The active control acts as an equality constraint, whereas the
* non-active controls should be interpreted as inequality
* constraints (upper or lower bounds). For instance, a PRODUCER's
* BHP constraint defines a minimum acceptable bottom-hole pressure
* value for the well.
/**
* Type of well control equation or inequality constraint.
*/
enum WellControlType {
BHP, /**< Well constrained by BHP target */
RESERVOIR_RATE, /**< Well constrained by reservoir volume flow rate */
SURFACE_RATE /**< Well constrained by surface volume flow rate */
};
/**
* Controls for a single well.
* Each control specifies a well rate or bottom-hole pressure. Only
* one control can be active at a time, indicated by current. The
* meaning of each control's target value depends on the control type:
*
* - BHP -> target pressure in Pascal.
* - RESERVOIR_RATE -> target reservoir volume rate in cubic(meter)/second
* - SURFACE_RATE -> target surface volume rate in cubic(meter)/second
*
* The sign convention for RATE targets is as follows:
*
* - (+) Fluid flowing into reservoir, i.e. injecting.
* - (-) Fluid flowing out of reservoir, i.e. producing.
*
* For *_RATE controls, the distribution of phases used for the control
* is also needed. For example, a total rate control should have 1.0
* for each phase, whereas a control on oil rate should have 1.0 for
* the oil phase and 0.0 for the rest. For BHP controls, this is unused.
* The active control acts as an equality constraint, whereas the
* non-active controls should be interpreted as inequality
* constraints (upper or lower bounds). For instance, a PRODUCER's
* BHP constraint defines a minimum acceptable bottom-hole pressure
* value for the well.
*/
struct WellControls
{
int num; /** Number of controls. */
enum WellControlType *type; /** Array of control types.*/
double *target; /** Array of control targets */
double *distr; /** Array of rate control distributions,
Wells::number_of_phases numbers per control */
int current; /** Index of current active control. */
/**
* Number of controls.
*/
int num;
void *data; /** Internal management structure. */
/**
* Array of control types.
*/
enum WellControlType *type;
/**
* Array of control targets
*/
double *target;
/**
* Array of rate control distributions,
* <CODE>Wells::number_of_phases</CODE> numbers for each control
*/
double *distr;
/**
* Index of current active control.
*/
int current;
/**
* Internal management structure.
*/
void *data;
};
/** Data structure aggregating static information about all wells in a scenario. */
/**
* Data structure aggregating static information about all wells in a scenario.
*/
struct Wells
{
int number_of_wells; /** Number of wells. */
int number_of_phases; /** Number of phases. */
int number_of_wells; /**< Number of wells. */
int number_of_phases; /**< Number of phases. */
enum WellType *type; /** Array of well types. */
double *depth_ref; /** Array of well bhp reference depths. */
double *comp_frac; /** Component fractions for each well, size is (number_of_wells*number_of_phases).
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
*/
int *well_connpos; /** Array of indices into well_cells (and WI).
* For a well w, well_connpos[w] and well_connpos[w+1] yield
* start and one-beyond-end indices into the well_cells array
* for accessing w's perforation cell indices.
*/
int *well_cells; /** Array of perforation cell indices.
* Size is number of perforations (== well_connpos[number_of_wells]).
*/
double *WI; /** Well productivity index, same size and structure as well_cells. */
struct WellControls **ctrls; /** Well controls, one set of controls for each well. */
/**
* Array of well types.
*/
enum WellType *type;
char **name; /** Well names. One string for each well. */
/**
* Array of well types.
*/
double *depth_ref;
void *data; /** Internal management structure. */
/**
* Component fractions for each well. Array of size
* <CODE>number_of_wells * number_of_phases</CODE>.
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
*/
double *comp_frac;
/**
* Array of indices into well_cells (and WI). For a well @c w,
* <CODE>well_connpos[w]</CODE> and <CODE>well_connpos[w+1]</CODE> are start
* and one-beyond-end indices into the @c well_cells array for accessing
* @c w's perforation cell indices.
*/
int *well_connpos;
/**
* Array of perforation cell indices.
* Size is number of perforations (== well_connpos[number_of_wells]).
*/
int *well_cells;
/**
* Well productivity index, same size and structure as well_cells.
*/
double *WI;
/**
* Well controls, one set of controls for each well.
*/
struct WellControls **ctrls;
/**
* Well names. One string for each well.
*/
char **name;
/**
* Internal management structure.
*/
void *data;
};
/** Data structure aggregating dynamic information about all wells in a scenario.
* All arrays in this structure contain data for each perforation,
* ordered the same as Wells::well_cells and Wells:WI. The array
* sizes are, respectively,
/**
* Data structure aggregating dynamic information about all wells in a scenario.
* All arrays in this structure contain data for each perforation, ordered the
* same as Wells::well_cells and Wells:WI. The array sizes are, respectively,
*
* gpot n*NP
* wdp NP
* A n²*NP (matrix in column-major (i.e., Fortran) order).
* phasemob n*NP
*
* in which "n" denotes the number of active fluid phases (and
* constituent components) and "NP" is the total number of
* perforations, <CODE>well_connpos[ number_of_wells ]</CODE>.
* in which "n" denotes the number of active fluid phases (and constituent
* components) and "NP" is the total number of perforations,
* <CODE>well_connpos[ number_of_wells ]</CODE>.
*/
struct CompletionData
{
double *gpot; /** Gravity potentials. */
double *A; /** Volumes to surface-components matrix, A = RB^{-1}. */
double *phasemob; /** Phase mobilities. */
/**
* Gravity potentials.
*/
double *wdp;
/**
* Volumes to surface-components matrix, A = RB^{-1}.
*/
double *A;
/**
* Phase mobilities for all perforations, stored consecutively with the
* phase index cycling the most rapidly.
*/
double *phasemob;
};
/**

32
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -216,9 +216,9 @@ namespace Opm
const int nperf = wells_->well_connpos[nw];
const int dim = grid_.dimensions;
const double grav = gravity_ ? gravity_[dim - 1] : 0.0;
wellperf_gpot_.clear();
wellperf_gpot_.resize(np*nperf, 0.0);
if (grav == 0.0) {
wellperf_wdp_.clear();
wellperf_wdp_.resize(nperf, 0.0);
if (not (std::abs(grav) > 0.0)) {
return;
}
@ -228,9 +228,9 @@ namespace Opm
// Main loop, iterate over all perforations,
// using the following formula (by phase):
// gpot(perf) = g*(perf_z - well_ref_z)*rho(perf)
// where the phase densities rho(perf) are taken to be
// the densities in the perforation cell.
// wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
// where the total density rho(perf) is taken to be
// sum_p (rho_p*saturation_p) in the perforation cell.
for (int w = 0; w < nw; ++w) {
const double ref_depth = wells_->depth_ref[w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w + 1]; ++j) {
@ -239,7 +239,8 @@ namespace Opm
props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &rho[0]);
for (int phase = 0; phase < np; ++phase) {
wellperf_gpot_[np*j + phase] = rho[phase]*grav*(cell_depth - ref_depth);
const double s_phase = state.saturation()[np*cell + phase];
wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
}
}
}
@ -478,10 +479,7 @@ namespace Opm
std::copy(cM, cM + np, wpM);
} else {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
double perf_p = bhp + wellperf_wdp_[j];
// Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to
@ -512,7 +510,7 @@ namespace Opm
const double* z = &state.surfacevol()[0];
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
completion_data.wdp = ! wellperf_wdp_.empty() ? &wellperf_wdp_[0] : 0;
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
cfs_tpfa_res_wells wells_tmp;
@ -599,7 +597,7 @@ namespace Opm
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
completion_data.wdp = ! wellperf_wdp_.empty() ? const_cast<double*>(&wellperf_wdp_[0]) : 0;
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
cfs_tpfa_res_wells wells_tmp;
@ -636,16 +634,10 @@ namespace Opm
// Compute well perforation pressures (not done by the C code).
if (wells_ != 0) {
const int nw = wells_->number_of_wells;
const int np = props_.numPhases();
for (int w = 0; w < nw; ++w) {
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
well_state.perfPress()[j] = perf_p;
well_state.perfPress()[j] = bhp + wellperf_wdp_[j];
}
}
}

2
opm/core/pressure/CompressibleTpfa.hpp
View File

@ -134,7 +134,7 @@ namespace Opm
struct cfs_tpfa_res_data* h_;
// ------ Data that will be modified for every solve. ------
std::vector<double> wellperf_gpot_;
std::vector<double> wellperf_wdp_;
std::vector<double> initial_porevol_;
// ------ Data that will be modified for every solver iteration. ------

32
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -420,15 +420,15 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
const double *wpress,
struct cfs_tpfa_res_impl *pimpl )
{
int c, i, w, np2;
double pw, dp;
double *WI, *gpot, *Ap, *pmobp;
double *pflux, *dpflux;
int c, i, w, np2;
double pw, dp;
const double *WI, *wdp, *Ap, *pmobp;
double *pflux, *dpflux, gpot[3] = { 0.0 };
struct Wells *W;
assert (wells->W != NULL);
assert (wells->W->number_of_phases <= 3);
W = wells->W;
@ -436,7 +436,7 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
assert (W->data != NULL);
WI = W->WI;
gpot = wells->data->gpot;
wdp = wells->data->wdp;
Ap = wells->data->A;
pmobp = wells->data->phasemob;
@ -449,11 +449,11 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
pw = wpress[w];
for (; i < W->well_connpos[w + 1]; i++,
gpot += np, Ap += np2, pmobp += np,
Ap += np2, pmobp += np,
pflux += np, dpflux += 2 * np) {
c = W->well_cells[i];
dp = pw - cpress[c];
dp = pw + wdp[i]- cpress[c];
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
pimpl->flux_work,
@ -909,8 +909,8 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
{
int w, i, c, np, np2, nc;
int is_neumann, is_open;
double pw, dp;
double *WI, *gpot, *pmobp;
double pw, dp, gpot[3] = { 0.0 };
const double *WI, *wdp, *pmobp;
const double *Ac, *dAc;
struct Wells *W;
@ -923,7 +923,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
W = wells->W;
WI = W->WI;
gpot = wells->data->gpot;
wdp = wells->data->wdp;
pmobp = wells->data->phasemob;
is_neumann = 1;
@ -932,14 +932,13 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
pw = wpress[ w ];
is_open = W->ctrls[w]->current >= 0;
for (; i < W->well_connpos[w + 1];
i++, gpot += np, pmobp += np) {
for (; i < W->well_connpos[w + 1]; i++, pmobp += np) {
c = W->well_cells[ i ];
Ac = cq->Ac + (c * np2);
dAc = cq->dAc + (c * np2);
dp = pw - cpress[ c ];
dp = pw + wdp[i] - cpress[ c ];
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
@ -1091,7 +1090,7 @@ compute_wflux(int np ,
{
int w, c, i, p;
double pw, dp, t;
const double *pmob;
const double *pmob, *wdp;
struct Wells *W;
struct CompletionData *cdata;
@ -1102,13 +1101,14 @@ compute_wflux(int np ,
W = wells->W;
cdata = wells->data;
wdp = cdata->wdp;
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[w];
for (; i < W->well_connpos[w + 1]; i++) {
c = W->well_cells[ i ];
dp = pw - cpress[c];
dp = pw + wdp[ i ] - cpress[c];
if (dp > 0) { pmob = cdata->phasemob + (i * np); } /* w->c */
else { pmob = pmobc + (c * np); } /* c->w */

238
opm/core/pressure/tpfa/cfs_tpfa_residual.h
View File

@ -25,6 +25,37 @@
#include <opm/core/pressure/tpfa/compr_source.h>
/**
* \file
* Public interface to assembler for (compressible) discrete pressure system
* based on two-point flux approximation method. The assembler implements a
* residual formulation that for a single cell \f$i\f$ reads
* \f[
* \mathit{pv}_i\cdot (1 - \sum_\alpha A_i^{-1}(p^{n+1},z^n)z^n) +
* \Delta t\sum_\alpha A_i^{-1}(p^{n+1},z^n) \Big(\sum_j A_{ij}v_{ij}^{n+1} -
* q_i^{n+1}\Big) = 0
* \f]
* in which \f$\mathit{pv}_i\f$ is the (constant or pressure-dependent)
* pore-volume of cell \f$i\f$. Moreover, \f$\Delta t\f$ is the time step size,
* \f$n\f$ denotes the time level, \f$A\f$ is the pressure and mass-dependent
* fluid matrix that converts phase volumes at reservoir conditions into
* component volumes at surface conditions and \f$v_{ij}\f$ is the vector of
* outward (with respect to cell \f$i\f$) phase fluxes across the \f$ij\f$
* cell-interface.
*
* This module's usage model is intended to be
* -# Construct assembler
* -# for (each time step)
* -# while (pressure not converged)
* -# Assemble appropriate (Jacobian) system of linear equations
* -# Solve Jacobian system to derive pressure increments
* -# Include increments into current state
* -# Check convergence
* -# Derive fluxes and, optionally, interface pressures
* -# Solve transport by some means
* -# Destroy assembler
*/
#ifdef __cplusplus
extern "C" {
#endif
@ -33,36 +64,121 @@ struct cfs_tpfa_res_impl;
struct CSRMatrix;
struct compr_quantities_gen;
/**
* Type encapsulating well topology and completion data (e.g., phase mobilities
* per connection (perforation)).
*/
struct cfs_tpfa_res_wells {
/**
* All wells pertaining to a particular linear system assembly.
* Must include current controls/targets and complete well topology.
*/
struct Wells *W ;
/**
* Completion data describing the fluid state at the current time level.
*/
struct CompletionData *data;
};
/**
* Type encapsulating all driving forces affecting the discrete pressure system.
*/
struct cfs_tpfa_res_forces {
struct cfs_tpfa_res_wells *wells;
struct compr_src *src ;
struct cfs_tpfa_res_wells *wells; /**< Wells */
struct compr_src *src ; /**< Explicit source terms */
};
/**
* Result structure that presents the fully assembled system of linear
* equations, linearised around the current pressure point.
*/
struct cfs_tpfa_res_data {
struct CSRMatrix *J;
double *F;
struct CSRMatrix *J; /**< Jacobian matrix */
double *F; /**< Residual vector (right-hand side) */
struct cfs_tpfa_res_impl *pimpl;
struct cfs_tpfa_res_impl *pimpl; /**< Internal management structure */
};
/**
* Construct assembler for system of linear equations.
*
* @param[in] G Grid
* @param[in] wells Well description. @c NULL in case of no wells.
* For backwards compatibility, the constructor also
* interprets <CODE>(wells != NULL) &&
* (wells->W == NULL)</CODE> as "no wells present", but new
* code should use <CODE>wells == NULL</CODE> to signify
* "no wells".
* @param[in] nphases Number of active fluid phases in this simulation run.
* Needed to correctly size various internal work arrays.
* @return Fully formed assembler structure suitable for forming systems of
* linear equations using, e.g., function cfs_tpfa_res_assemble(). @c NULL in
* case of allocation failure. Must be destroyed using function
* cfs_tpfa_res_destroy().
*/
struct cfs_tpfa_res_data *
cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_wells *wells ,
int nphases);
/**
* Destroy assembler for system of linear equations.
*
* Disposes of all resources acquired in a previous call to construction
* function cfs_tpfa_res_construct(). Note that the statement
* <CODE>
* cfs_tpfa_res_destroy(NULL)
* </CODE>
* is supported and benign (i.e., behaves like <CODE>free(NULL)</CODE>).
*
* @param[in,out] h On input - assembler obtained through a previous call to
* construction function cfs_tpfa_res_construct(). On output - invalid pointer.
*/
void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible and
there are no pressure conditions on wells or boundaries.
Otherwise return 0. */
/**
* Assemble system of linear equations by linearising the residual around the
* current pressure point. Assume incompressible rock (i.e., that the
* pore-volume is independent of pressure).
*
* The fully assembled system is presented in <CODE>h->J</CODE> and
* <CODE>h->F</CODE> and must be solved separately using external software.
*
* @param[in] G Grid.
* @param[in] dt Time step size \f$\Delta t\f$.
* @param[in] forces Driving forces.
* @param[in] zc Component volumes, per pore-volume, at surface
* conditions for all components in all cells stored
* consecutively per cell. Array of size
* <CODE>G->number_of_cells * cq->nphases</CODE>.
* @param[in] cq Compressible quantities describing the current fluid
* state. Fields @c Ac, @c dAc, @c Af, and
* @c phasemobf must be valid.
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute().
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressures. One scalar value per grid cell.
* @param[in] wpress Well (bottom-hole) pressures. One scalar value per
* well. @c NULL in case of no wells.
* @param[in] porevol Pore-volumes. One (positive) scalar value for each
* grid cell.
* @param[in,out] h On input-a valid (non-@c NULL) assembler obtained
* from a previous call to constructor function
* cfs_tpfa_res_construct(). On output-valid assembler
* that includes the new system of linear equations in
* its @c J and @c F fields.
*
* @return 1 if the assembled matrix was adjusted to remove a singularity. This
* happens if all fluids are incompressible and there are no pressure conditions
* on wells or boundaries. Otherwise return 0.
*/
int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
double dt,
@ -76,11 +192,48 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
const double *porevol,
struct cfs_tpfa_res_data *h);
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible, the
rock is incompressible, and there are no pressure conditions on
wells or boundaries.
Otherwise return 0. */
/**
* Assemble system of linear equations by linearising the residual around the
* current pressure point. Assume compressible rock (i.e., that the pore-volume
* depends on pressure).
*
* The fully assembled system is presented in <CODE>h->J</CODE> and
* <CODE>h->F</CODE> and must be solved separately using external software.
*
* @param[in] G Grid.
* @param[in] dt Time step size \f$\Delta t\f$.
* @param[in] forces Driving forces.
* @param[in] zc Component volumes, per pore-volume, at surface
* conditions for all components in all cells stored
* consecutively per cell. Array of size
* <CODE>G->number_of_cells * cq->nphases</CODE>.
* @param[in] cq Compressible quantities describing the current fluid
* state. Fields @c Ac, @c dAc, @c Af, and
* @c phasemobf must be valid.
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute().
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressures. One scalar value per grid cell.
* @param[in] wpress Well (bottom-hole) pressures. One scalar value per
* well. @c NULL in case of no wells.
* @param[in] porevol Pore-volumes. One (positive) scalar value for each
* grid cell.
* @param[in] porevol0 Pore-volumes at start of time step (i.e., at time
* level \f$n\f$). One (positive) scalar value for
* each grid cell.
* @param[in] rock_comp Rock compressibility. One non-negative scalar for
* each grid cell.
* @param[in,out] h On input-a valid (non-@c NULL) assembler obtained
* from a previous call to constructor function
* cfs_tpfa_res_construct(). On output-valid assembler
* that includes the new system of linear equations in
* its @c J and @c F fields.
*
* @return 1 if the assembled matrix was adjusted to remove a singularity. This
* happens if all fluids are incompressible, the rock is incompressible, and
* there are no pressure conditions on wells or boundaries. Otherwise return 0.
*/
int
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G,
@ -97,6 +250,38 @@ cfs_tpfa_res_comprock_assemble(
const double *rock_comp,
struct cfs_tpfa_res_data *h);
/**
* Derive interface (total) Darcy fluxes from (converged) pressure solution.
*
* @param[in] G Grid
* @param[in] forces Driving forces. Must correspond to those used when
* forming the system of linear equations, e.g., in the
* call to function cfs_tpfa_res_assemble().
* @param[in] np Number of fluid phases (and components).
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute(). Must correspond to equally
* named parameter of the assembly functions.
* @param[in] pmobc Phase mobilities stored consecutively per cell with
* phase index cycling the most rapidly. Array of size
* <CODE>G->number_of_cells * np</CODE>.
* @param[in] pmobf Phase mobilities stored consecutively per interface
* with phase index cycling the most rapidly. Array of
* size <CODE>G->number_of_faces * np</CODE>.
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressure values. One (positive) scalar for each
* grid cell.
* @param[in] wpress Well (bottom-hole) pressure values. One (positive)
* scalar value for each well. @c NULL in case of no
* wells.
* @param[out] fflux Total Darcy interface fluxes. One scalar value for
* each interface (inter-cell connection). Array of size
* <CODE>G->number_of_faces</CODE>.
* @param[out] wflux Total Darcy well connection fluxes. One scalar value
* for each well connection (perforation). Array of size
* <CODE>forces->wells->W->well_connpos
* [forces->wells->W->number_of_wells]</CODE>.
*/
void
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_forces *forces ,
@ -110,6 +295,31 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
double *fflux ,
double *wflux );
/**
* Derive interface pressures from (converged) pressure solution.
*
* @param[in] G Grid
* @param[in] np Number of fluid phases (and components).
* @param[in] htrans Background one-sided ("half") transmissibilities as
* defined by function tpfa_htrans_compute().
* @param[in] pmobf Phase mobilities stored consecutively per interface
* with phase index cycling the most rapidly. Array of
* size <CODE>G->number_of_faces * np</CODE>.
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] h System assembler. Must correspond to the assembler
* state used to form the final system of linear equations
* from which the converged pressure solution was derived.
* @param[in] cpress Cell pressure values. One (positive) scalar for each
* grid cell.
* @param[in] fflux Total Darcy interface fluxes. One scalar value for
* each interface (inter-cell connection). Array of size
* <CODE>G->number_of_faces</CODE>. Typically computed
* using function cfs_tpfa_res_flux().
* @param[out] fpress Interface pressure values. One (positive) scalar for
* each interface. Array of size
* <CODE>G->number_of_faces</CODE>.
*/
void
cfs_tpfa_res_fpress(struct UnstructuredGrid *G,
int np,

98
opm/core/pressure/tpfa/compr_quant_general.h
View File

@ -22,23 +22,109 @@
#include <stddef.h>
/**
* \file
* Module defining derived fluid quantities needed to discretise compressible
* and miscible pressure (flow) problems.
*/
#ifdef __cplusplus
extern "C" {
#endif
/**
* Aggregate structure that represents an atomic view of the current fluid
* state. These quantities are used directly in the cfs_tpfa_residual module to
* discretise a particular, linearised flow problem.
*/
struct compr_quantities_gen {
int nphases; /* Number of phases/components */
/**
* Number of fluid phases. The pressure solvers also assume that the number
* of fluid components (at surface conditions) equals the number of fluid
* phases.
*/
int nphases;
double *Ac; /* RB^{-1} per cell */
double *dAc; /* d/dp (RB^{-1}) per cell */
double *Af; /* RB^{-1} per face */
double *phasemobf; /* Phase mobility per face */
double *voldiscr; /* Volume discrepancy per cell */
/**
* Pressure and mass-dependent fluid matrix that converts phase volumes at
* reservoir conditions into component volumes at surface conditions. Obeys
* the defining formula
* \f[
* A = RB^{-1}
* \f]
* in which \f$R\f$ denotes the miscibility ratios (i.e., the dissolved
* gas-oil ratio, \f$R_s\f$ and the evaporated oil-gas ratio, \f$R_v\f$)
* collected as a matrix and \f$B\f$ is the diagonal matrix of
* formation-volume expansion factors. The function is sampled in each grid
* cell. Array of size <CODE>nphases * nphases * nc</CODE>.
*/
double *Ac;
/**
* Derivative of \f$A\f$ with respect to pressure,
* \f[
* \frac{\partial A}{\partial p} = \frac{\partial R}{\partial p}B^{-1} +
* R\frac{\partial B^{-1}}{\partial p} = (\frac{\partial R}{\partial p} -
* A\frac{\partial B}{\partial p})B^{-1}
* \f]
* sampled in each grid cell. Array of size
* <CODE>nphases * nphases * nc</CODE>.
*/
double *dAc;
/**
* Fluid matrix sampled at each interface. Possibly as a result of an
* upwind calculation. Array of size <CODE>nphases * nphases * nf</CODE>.
*/
double *Af;
/**
* Phase mobility per interface. Possibly defined through an upwind
* calculation. Array of size <CODE>nphases * nf</CODE>.
*/
double *phasemobf;
/**
* Deceptively named "volume-discrepancy" term. Array of size @c nc.
* Unused by the cfs_tpfa_residual module.
*/
double *voldiscr;
};
/**
* Create management structure capable of storing derived fluid quantities
* pertaining to a reservoir model of a particular size.
*
* The resources acquired using function compr_quantities_gen_allocate() must
* be released using the destructor function compr_quantities_gen_deallocate().
*
* @param[in] nc Number of grid cells
* @param[in] nf Number of grid interfaces
* @param[in] np Number of fluid phases (and components)
*
* @return Fully formed management structure. Returns @c NULL in case of
* allocation failure.
*/
struct compr_quantities_gen *
compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
/**
* Release resources acquired in a previous call to constructor function
* compr_quantities_gen_allocate().
*
* Note that
* <CODE>
* compr_quantities_gen_deallocate(NULL)
* </CODE>
* is supported and benign (i.e., this statement behaves as
* <CODE>free(NULL)</CODE>).
*
* @param[in,out] cq On input - compressible quantity management structure
* obtained through a previous call to construction function
* compr_quantities_gen_allocate(). On output - invalid pointer.
*/
void
compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);

101
opm/core/pressure/tpfa/compr_source.h
View File

@ -36,28 +36,112 @@
#ifndef OPM_COMPR_SOURCE_H_HEADER
#define OPM_COMPR_SOURCE_H_HEADER
/**
* \file
* Data structures and support routines needed to represent explicit,
* compressible source terms.
*/
#ifdef __cplusplus
extern "C" {
#endif
/**
* Collection of explicit, compressible source terms.
*/
struct compr_src {
int nsrc;
int cpty;
/**
* Number of source terms.
*/
int nsrc;
int nphases;
/**
* Source term capacity. Client code should treat this member as read-only.
* The field is used in internal memory management.
*/
int cpty;
int *cell;
/**
* Number of fluid phases.
*/
int nphases;
/**
* Cells influenced by explicit source terms. Array of size @c cpty, the
* @c nsrc first elements (only) of which are valid.
*/
int *cell;
/**
* Total Darcy rate of inflow (measured at reservoir conditions) of each
* individual source term. Sign convention: Positive rate into reservoir
* (i.e., injection) and negative rate out of reservoir (production).
* Array of size @c cpty, the @c nsrc first elements (only) of which are
* valid.
*/
double *flux;
/**
* Injection composition for all explicit source terms. Not referenced for
* production sources (i.e., those terms for which <CODE>->flux[]</CODE> is
* negative). Array of size <CODE>nphases * cpty</CODE>, the
* <CODE>nphases * nsrc</CODE> of which (only) are valid.
*/
double *saturation;
};
/**
* Create a management structure that is, initially, capable of storing a
* specified number of sources defined by a particular number a fluid phases.
*
* @param[in] np Number of fluid phases. Must be positive to actually
* allocate any sources.
* @param[in] nsrc Initial management capacity. If positive, attempt to
* allocate that number of source terms. Otherwise, the
* initial capacity is treated as (and set to) zero.
* @return Fully formed management structure if <CODE>np > 0</CODE> and
* allocation success. @c NULL otherwise. The resources must be released using
* destructor function compr_src_deallocate().
*/
struct compr_src *
compr_src_allocate(int np, int nsrc);
/**
* Release memory resources acquired in a previous call to constructor function
* compr_src_allocate() and, possibly, source term insertion function
* append_compr_source_term().
*
* @param[in,out] src On input - source term management structure obtained
* through a previous call to construction function compr_src_allocate().
* On output - invalid pointer.
*/
void
compr_src_deallocate(struct compr_src *src);
/**
* Insert a new explicit source term into an existing collection.
*
* @param[in] c Cell influenced by this source term.
* @param[in] np Number of fluid phases. Used for consistency checking
* only.
* @param[in] v Source term total reservoir volume Darcy flux. Positive
* if the source term is to be interpreted as an injection
* source and negative otherwise.
* @param[in] sat Injection composition for this source term. Array of size
* @c np. Copied to internal storage, but the actual numbers
* are not inspected unless <CODE>v > 0.0</CODE>.
* @param[in,out] src On input - source term management structure obtained
* through a previous call to construction function compr_src_allocate() and,
* possibly, another call to function append_compr_source_term(). On output -
* source term collection that includes the new source term if successful and
* unchanged if not.
*
* @return One (1, true) if successful (i.e., if the new source term could be
* included in the existing collection) and zero (0, false) if not.
*/
int
append_compr_source_term(int c ,
int np ,
@ -65,6 +149,15 @@ append_compr_source_term(int c ,
const double *sat,
struct compr_src *src);
/**
* Empty source term collection while maintaining existing capacity.
*
* @param[in,out] src On input - an existing source term collection with a
* given number of sources and source capacity. On output - an empty source
* term collection (i.e., <CODE>src->nsrc == 0</CODE>) with an unchanged
* capacity.
*/
void
clear_compr_source_term(struct compr_src *src);