diff --git a/opm/core/newwells.h b/opm/core/newwells.h index 20cfce61..331d554e 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -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, + * Wells::number_of_phases 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 + * 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. + */ + double *comp_frac; + + /** + * Array of indices into well_cells (and WI). For a well @c w, + * well_connpos[w] and well_connpos[w+1] 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, well_connpos[ number_of_wells ]. + * in which "n" denotes the number of active fluid phases (and constituent + * components) and "NP" is the total number of perforations, + * well_connpos[ number_of_wells ]. */ 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; }; /** diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 8edb41a0..742a8a27 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -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(&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(&grid_); CompletionData completion_data; - completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast(&wellperf_gpot_[0]) : 0; + completion_data.wdp = ! wellperf_wdp_.empty() ? const_cast(&wellperf_wdp_[0]) : 0; completion_data.A = ! wellperf_A_.empty() ? const_cast(&wellperf_A_[0]) : 0; completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast(&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]; } } } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 0b54178e..bf9e2901 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -134,7 +134,7 @@ namespace Opm struct cfs_tpfa_res_data* h_; // ------ Data that will be modified for every solve. ------ - std::vector wellperf_gpot_; + std::vector wellperf_wdp_; std::vector initial_porevol_; // ------ Data that will be modified for every solver iteration. ------ diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index bcbacec2..cb3d8173 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -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 */ diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 5eddeed6..368cce10 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -25,6 +25,37 @@ #include +/** + * \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 (wells != NULL) && + * (wells->W == NULL) as "no wells present", but new + * code should use wells == NULL 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 + * + * cfs_tpfa_res_destroy(NULL) + * + * is supported and benign (i.e., behaves like free(NULL)). + * + * @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 h->J and + * h->F 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 + * G->number_of_cells * cq->nphases. + * @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 h->J and + * h->F 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 + * G->number_of_cells * cq->nphases. + * @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 + * G->number_of_cells * np. + * @param[in] pmobf Phase mobilities stored consecutively per interface + * with phase index cycling the most rapidly. Array of + * size G->number_of_faces * np. + * @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 + * G->number_of_faces. + * @param[out] wflux Total Darcy well connection fluxes. One scalar value + * for each well connection (perforation). Array of size + * forces->wells->W->well_connpos + * [forces->wells->W->number_of_wells]. + */ 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 G->number_of_faces * np. + * @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 + * G->number_of_faces. Typically computed + * using function cfs_tpfa_res_flux(). + * @param[out] fpress Interface pressure values. One (positive) scalar for + * each interface. Array of size + * G->number_of_faces. + */ void cfs_tpfa_res_fpress(struct UnstructuredGrid *G, int np, diff --git a/opm/core/pressure/tpfa/compr_quant_general.h b/opm/core/pressure/tpfa/compr_quant_general.h index 9fd13944..7027ee39 100644 --- a/opm/core/pressure/tpfa/compr_quant_general.h +++ b/opm/core/pressure/tpfa/compr_quant_general.h @@ -22,23 +22,109 @@ #include +/** + * \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 nphases * nphases * nc. + */ + 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 + * nphases * nphases * nc. + */ + double *dAc; + + /** + * Fluid matrix sampled at each interface. Possibly as a result of an + * upwind calculation. Array of size nphases * nphases * nf. + */ + double *Af; + + /** + * Phase mobility per interface. Possibly defined through an upwind + * calculation. Array of size nphases * nf. + */ + 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 + * + * compr_quantities_gen_deallocate(NULL) + * + * is supported and benign (i.e., this statement behaves as + * free(NULL)). + * + * @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); diff --git a/opm/core/pressure/tpfa/compr_source.h b/opm/core/pressure/tpfa/compr_source.h index d88c1b55..9beb0964 100644 --- a/opm/core/pressure/tpfa/compr_source.h +++ b/opm/core/pressure/tpfa/compr_source.h @@ -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 ->flux[] is + * negative). Array of size nphases * cpty, the + * nphases * nsrc 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 np > 0 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 v > 0.0. + * @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., src->nsrc == 0) with an unchanged + * capacity. + */ void clear_compr_source_term(struct compr_src *src);