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

Switch to a sparse/compressed boundary condition representation.

Browse Source 
Specifically, replace the existing flowbc_t (that was densely
represented on each interface, including internal interface and
external no-flow interfaces) with a new structure given by

    struct FlowBoundaryConditions

The semantics of this structure mirror those of "struct Wells" from
<opm/core/newwells.h>, but is currently mostly intended for simple,
incompressible flow purposes.

Update pressure solvers supporting boundary conditions to accommodate
the new boundary condition representation in the process.
This commit is contained in:
Bård Skaflestad 2012-03-06 20:07:35 +01:00
parent 7149d10cb2
commit 69bc8e16d6
9 changed files with 409 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
opm/core/pressure/cfsh.c
View File

@ -34,7 +34,7 @@
/* ---------------------------------------------------------------------- */
static int
cfsh_assemble_grid(flowbc_t *bc,
cfsh_assemble_grid(struct FlowBoundaryConditions *bc,
const double *Binv,
const double *gpress,
const double *src,
@ -126,8 +126,8 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W)
ngconn_tot = G->cell_facepos[nc];
fsh_define_linsys_arrays(new);
fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn,
W, new->pimpl);
fsh_define_impl_arrays(nc, G->number_of_faces, nnu,
ngconn_tot, new->max_ngconn, W, new->pimpl);
new->pimpl->sys = hybsys_allocate_unsymm(new->max_ngconn,
nc, ngconn_tot);
@ -177,7 +177,7 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W)
*/
/* ---------------------------------------------------------------------- */
void
cfsh_assemble(flowbc_t *bc,
cfsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *Biv,

174
opm/core/pressure/flow_bc.c
View File

@ -1,5 +1,5 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
@ -17,54 +17,172 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdlib.h>
#include "flow_bc.h"
#include <opm/core/pressure/flow_bc.h>
/* Create structure to hold flow boundary conditions for 'nf' faces.
*
* Return fully allocated structure or NULL in case of allocation
* failure. */
/* ---------------------------------------------------------------------- */
flowbc_t *
allocate_flowbc(size_t nf)
/* Compute an appropriate array dimension to minimise total number of
* (re-)allocations. */
/* ---------------------------------------------------------------------- */
static size_t
alloc_size(size_t n, size_t c)
/* ---------------------------------------------------------------------- */
{
size_t i;
flowbc_t *new;
if (c < n) {
c *= 2; /* log_2(n) allocations */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->type = malloc(nf * sizeof *new->type);
new->bcval = malloc(nf * sizeof *new->bcval);
if ((new->type == NULL) || (new->bcval == NULL)) {
deallocate_flowbc(new);
new = NULL;
} else {
for (i = 0; i < nf; i++) {
new->type [i] = UNSET;
new->bcval[i] = 0.0;
}
if (c < n) {
c = n; /* Typically for the first few allocs */
}
}
return new;
return c;
}
/* Release memory resources for dynamically allocated flow boundary
* condition structure. */
/* ---------------------------------------------------------------------- */
/* Put structure in a well-defined, initial state */
/* ---------------------------------------------------------------------- */
static void
initialise_structure(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
fbc->nbc = 0;
fbc->cpty = 0;
fbc->type = NULL;
fbc->value = NULL;
fbc->face = NULL;
}
/* ---------------------------------------------------------------------- */
static int
expand_tables(size_t nbc,
struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
int ok;
size_t alloc_sz;
void *p1, *p2, *p3;
ok = nbc <= fbc->cpty;
if (! ok) {
alloc_sz = alloc_size(nbc, fbc->cpty);
p1 = realloc(fbc->type , alloc_sz * sizeof *fbc->type );
p2 = realloc(fbc->value, alloc_sz * sizeof *fbc->value);
p3 = realloc(fbc->face , alloc_sz * sizeof *fbc->face );
ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL);
if (ok) {
fbc->type = p1;
fbc->value = p2;
fbc->face = p3;
fbc->cpty = alloc_sz;
} else {
free(p3); free(p2); free(p1);
}
}
return ok;
}
/* ======================================================================
* Public interface below separator
* ====================================================================== */
/* ---------------------------------------------------------------------- */
/* Allocate a 'FlowBoundaryConditions' structure, initially capable of
* managing 'nbc' individual boundary conditions. */
/* ---------------------------------------------------------------------- */
struct FlowBoundaryConditions *
flow_conditions_construct(size_t nbc)
/* ---------------------------------------------------------------------- */
{
int ok;
struct FlowBoundaryConditions *fbc;
fbc = malloc(1 * sizeof *fbc);
if (fbc != NULL) {
initialise_structure(fbc);
ok = expand_tables(nbc, fbc);
if (! ok) {
flow_conditions_destroy(fbc);
fbc = NULL;
}
}
return fbc;
}
/* ---------------------------------------------------------------------- */
/* Release memory resources managed by 'fbc', including the containing
* 'struct' pointer, 'fbc'. */
/* ---------------------------------------------------------------------- */
void
deallocate_flowbc(flowbc_t *fbc)
flow_conditions_destroy(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
if (fbc != NULL) {
free(fbc->bcval);
free(fbc->face );
free(fbc->value);
free(fbc->type );
}
free(fbc);
}
/* ---------------------------------------------------------------------- */
/* Append a new boundary condition to existing set.
*
* Return one (1) if successful, and zero (0) otherwise. */
/* ---------------------------------------------------------------------- */
int
flow_conditions_append(enum FlowBCType type ,
int face ,
double value,
struct FlowBoundaryConditions *fbc )
/* ---------------------------------------------------------------------- */
{
int ok;
ok = expand_tables(fbc->nbc + 1, fbc);
if (ok) {
fbc->type [ fbc->nbc ] = type ;
fbc->value[ fbc->nbc ] = value;
fbc->face [ fbc->nbc ] = face ;
fbc->nbc += 1;
}
return ok;
}
/* ---------------------------------------------------------------------- */
/* Clear existing set of boundary conditions */
/* ---------------------------------------------------------------------- */
void
flow_conditions_clear(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
assert (fbc != NULL);
fbc->nbc = 0;
}

48
opm/core/pressure/flow_bc.h
View File

@ -22,26 +22,56 @@
#include <stddef.h>
#ifdef __cplusplus
extern "C" {
#endif
enum flowbc_type { UNSET, PRESSURE, FLUX };
enum FlowBCType { BC_NOFLOW ,
BC_PRESSURE ,
BC_FLUX_TOTVOL };
typedef struct {
enum flowbc_type *type;
double *bcval;
} flowbc_t;
/* Boundary condition structure.
*
* Condition i (in [0 .. nbc-1]) affects (outer) interface face[i], is
* of type type[i], and specifies a target value of value[i].
*
* The field 'cpty' is for internal use by the implementation. */
struct FlowBoundaryConditions {
size_t nbc; /* Current number of bdry. conditions */
size_t cpty; /* Internal management. Do not touch */
enum FlowBCType *type; /* Condition type */
double *value; /* Condition value (target) */
int *face; /* Outer faces affected by ind. target */
};
flowbc_t *
allocate_flowbc(size_t nf);
/* Allocate a 'FlowBoundaryConditions' structure, initially capable of
* managing 'nbc' individual boundary conditions. */
struct FlowBoundaryConditions *
flow_conditions_construct(size_t nbc);
/* Release memory resources managed by 'fbc', including the containing
* 'struct' pointer, 'fbc'. */
void
deallocate_flowbc(flowbc_t *fbc);
flow_conditions_destroy(struct FlowBoundaryConditions *fbc);
/* Append a new boundary condition to existing set.
*
* Return one (1) if successful, and zero (0) otherwise. */
int
flow_conditions_append(enum FlowBCType type ,
int face ,
double value,
struct FlowBoundaryConditions *fbc );
/* Clear existing set of boundary conditions */
void
flow_conditions_clear(struct FlowBoundaryConditions *fbc);
#ifdef __cplusplus
}
#endif

4
opm/core/pressure/fsh.h
View File

@ -79,7 +79,7 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W);
/** Assembles the hybridized linear system for face pressures.
*/
void
cfsh_assemble(flowbc_t *bc,
cfsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *Biv,
@ -132,7 +132,7 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W);
* be modified). Must already be constructed.
*/
void
ifsh_assemble(flowbc_t *bc,
ifsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *gpress,

49
opm/core/pressure/fsh_common_impl.c
View File

@ -157,23 +157,28 @@ fsh_count_grid_dof(struct UnstructuredGrid *G, int *max_ngdof, size_t *sum_ngdof
/* Impose boundary conditions on local contribution to global system. */
/* ---------------------------------------------------------------------- */
int
fsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
fsh_impose_bc(int nconn, int *conn, struct FlowBoundaryConditions *bc,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int i, npp, f;
int i, j, npp, f;
npp = 0;
for (i = 0; i < nconn; i++) {
f = conn[i];
if (bc->type[f] == PRESSURE) {
pimpl->work [npp] = bc->bcval[f];
pimpl->iwork[npp] = i;
j = pimpl->bdry_condition[ f ];
npp += 1;
} else if (bc->type[f] == FLUX) {
pimpl->sys->r[i] -= bc->bcval[f];
if (j != -1) {
if (bc->type[j] == BC_PRESSURE) {
pimpl->work [npp] = bc->value[j];
pimpl->iwork[npp] = i;
npp += 1;
} else if (bc->type[j] == BC_FLUX_TOTVOL) {
pimpl->sys->r[i] -= bc->value[j];
}
}
}
@ -189,9 +194,32 @@ fsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
}
/* ---------------------------------------------------------------------- */
void
fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int f, i;
for (i = 0; i < pimpl->nf; i++) { pimpl->bdry_condition[ i ] = -1; }
if (fbc != NULL) {
for (i = 0; ((size_t) i) < fbc->nbc; i++) {
f = fbc->face[ i ];
assert ((0 <= f) && (f < pimpl->nf));
pimpl->bdry_condition[ f ] = i;
}
}
}
/* ---------------------------------------------------------------------- */
void
fsh_define_impl_arrays(size_t nc,
size_t nf,
size_t nnu,
size_t nhf,
size_t max_ncf,
@ -206,8 +234,10 @@ fsh_define_impl_arrays(size_t nc,
pimpl->gdof = pimpl->gdof_pos + (nc + 1);
pimpl->iwork = pimpl->gdof + nhf;
pimpl->bdry_condition = pimpl->iwork + max_ncf;
if (W != NULL) {
pimpl->cwell_pos = pimpl->iwork + max_ncf;
pimpl->cwell_pos = pimpl->bdry_condition + nf;
pimpl->cwells = pimpl->cwell_pos + nc + 1;
pimpl->WI = pimpl->work + max_ncf;
@ -256,6 +286,7 @@ fsh_compute_table_sz(struct UnstructuredGrid *G, well_t *W, int max_ngconn,
*idata_sz = nc + 1; /* gdof_pos */
*idata_sz += ngconn_tot; /* gdof */
*idata_sz += G->number_of_faces; /* bdry_condition */
*idata_sz += max_ngconn; /* iwork */
*ddata_sz = 2 * (*nnu); /* rhs + soln */

9
opm/core/pressure/fsh_common_impl.h
View File

@ -44,6 +44,8 @@ struct fsh_impl {
double *work; /* Scratch array, floating point */
int *iwork; /* Scratch array, integers */
int *bdry_condition; /* Map face->boundary condition ID */
/* Linear storage goes here... */
int *idata; /* Actual storage array, integers */
double *ddata; /* Actual storage array, floating point */
@ -56,14 +58,19 @@ fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz);
void
fsh_count_grid_dof(struct UnstructuredGrid *G, int *max_ngdof, size_t *sum_ngdof2);
void
fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
struct fsh_impl *pimpl);
int
fsh_impose_bc(int ndof,
int *dof,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
struct fsh_impl *pimpl);
void
fsh_define_impl_arrays(size_t nc,
size_t nf,
size_t nnu,
size_t nhf,
size_t max_ncf,

28
opm/core/pressure/ifsh.c
View File

@ -60,7 +60,7 @@ ifsh_set_effective_well_params(const double *WI,
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_grid(flowbc_t *bc,
ifsh_assemble_grid(struct FlowBoundaryConditions *bc,
const double *Binv,
const double *gpress,
const double *src,
@ -100,12 +100,12 @@ ifsh_assemble_grid(flowbc_t *bc,
/* ---------------------------------------------------------------------- */
static void
ifsh_impose_well_control(int c,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int ngconn, nwconn, i, w1, w2, wg, f;
int ngconn, nwconn, i, j, w1, w2, wg, f;
int *pgconn, *gconn, *pwconn, *wconn;
double bhp;
@ -130,11 +130,12 @@ ifsh_impose_well_control(int c,
/* Adapt local system to prescribed boundary pressures (r->w) */
for (i = 0; i < ngconn; i++) {
f = gconn[i];
j = ifsh->pimpl->bdry_condition[ f ];
if (bc->type[f] == PRESSURE) {
if (j != -1) {
for (w1 = 0; w1 < nwconn; w1++) {
/* Eliminate prescribed (boundary) pressure value */
r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->bcval[f];
r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->value[j];
r2w[i + w1*ngconn] = 0.0;
}
@ -152,8 +153,13 @@ ifsh_impose_well_control(int c,
/* Well->reservoir */
for (i = 0; i < ngconn; i++) {
assert ((bc->type[gconn[i]] != PRESSURE) ||
#ifndef NDEBUG
j = ifsh->pimpl->bdry_condition[ gconn[i] ];
assert ((j == -1) ||
(bc->type[j] != BC_PRESSURE) ||
!(fabs(r2w[i + w1*ngconn]) > 0.0));
#endif
r [i] -= r2w[i + w1*ngconn] * bhp;
r2w[i + w1*ngconn] = 0.0;
@ -177,7 +183,7 @@ ifsh_impose_well_control(int c,
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_well(flowbc_t *bc,
ifsh_assemble_well(struct FlowBoundaryConditions *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
@ -291,8 +297,8 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W)
ngconn_tot = G->cell_facepos[nc];
fsh_define_linsys_arrays(new);
fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn,
W, new->pimpl);
fsh_define_impl_arrays(nc, G->number_of_faces, nnu,
ngconn_tot, new->max_ngconn, W, new->pimpl);
new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn,
nc, ngconn_tot);
@ -343,7 +349,7 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W)
* pressure gpress, boundary conditions bc, and source terms src. */
/* ---------------------------------------------------------------------- */
void
ifsh_assemble(flowbc_t *bc,
ifsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *gpress,
@ -355,6 +361,8 @@ ifsh_assemble(flowbc_t *bc,
{
int npp; /* Number of prescribed pressure values */
fsh_map_bdry_condition(bc, ifsh->pimpl);
hybsys_schur_comp_symm(ifsh->pimpl->nc,
ifsh->pimpl->gdof_pos,
Binv, ifsh->pimpl->sys);

199
opm/core/pressure/tpfa/cfs_tpfa.c
View File

@ -36,6 +36,9 @@ struct densrat_util {
struct cfs_tpfa_impl {
/* Mapping (cell,face) -> half-face */
int *f2hf;
/* Reservoir flow */
double *ctrans;
double *accum;
@ -140,7 +143,7 @@ static struct cfs_tpfa_impl *
impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
/* ---------------------------------------------------------------------- */
{
size_t nnu, ngconn, nwperf;
size_t nnu, nf, ngconn, nwperf;
struct cfs_tpfa_impl *new;
size_t ddata_sz;
@ -158,10 +161,11 @@ impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
ddata_sz = 2 * nnu; /* b, x */
/* Reservoir */
nf = G->number_of_faces;
ddata_sz += 1 * ngconn; /* ctrans */
ddata_sz += 1 * G->number_of_cells; /* accum */
ddata_sz += np * G->number_of_faces; /* masstrans_f */
ddata_sz += np * G->number_of_faces; /* gravtrans_f */
ddata_sz += np * nf; /* masstrans_f */
ddata_sz += np * nf; /* gravtrans_f */
/* Wells */
ddata_sz += 1 * nwperf; /* wtrans */
@ -169,15 +173,17 @@ impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
ddata_sz += np * nwperf; /* masstrans_p */
ddata_sz += np * nwperf; /* gravtrans_p */
ddata_sz += 1 * G->number_of_faces; /* scratch_f */
ddata_sz += 1 * nf; /* scratch_f */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->f2hf = malloc(2 * nf * sizeof *new->f2hf );
new->ddata = malloc(ddata_sz * sizeof *new->ddata);
new->ratio = allocate_densrat(G, W, np);
if (new->ddata == NULL || new->ratio == NULL) {
if (new->f2hf == NULL ||
new->ddata == NULL || new->ratio == NULL) {
impl_deallocate(new);
new = NULL;
}
@ -392,7 +398,7 @@ set_dynamic_trans(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
static void
set_dynamic_grav(struct UnstructuredGrid *G ,
flowbc_t *bc ,
struct FlowBoundaryConditions *bc ,
const double *trans ,
const double *gravcap_f,
struct compr_quantities *cq ,
@ -405,7 +411,7 @@ set_dynamic_grav(struct UnstructuredGrid *G ,
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if (((c1 >= 0) && (c2 >= 0)) || (bc->type[f] == PRESSURE)) {
if ((c1 >= 0) && (c2 >= 0)) {
for (p = 0; p < cq->nphases; p++, i++) {
ratio->x[i] = trans[f] * gravcap_f[i] * cq->phasemobf[i];
}
@ -415,6 +421,18 @@ set_dynamic_grav(struct UnstructuredGrid *G ,
}
}
}
if (bc != NULL) {
for (i = 0; ((size_t) i) < bc->nbc; i++) {
if (bc->type[ i ] == BC_PRESSURE) {
f = bc->face[ i ];
for (p = cq->nphases * f; p < cq->nphases * (f + 1); p++) {
ratio->x[p] = trans[f] * gravcap_f[p] * cq->phasemobf[p];
}
}
}
}
}
@ -561,7 +579,7 @@ static void
compute_psys_contrib(struct UnstructuredGrid *G,
well_t *W,
struct completion_data *wdata,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
struct compr_quantities *cq,
double dt,
const double *trans,
@ -626,20 +644,16 @@ compute_psys_contrib(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
static int
assemble_cell_contrib(struct UnstructuredGrid *G,
flowbc_t *bc,
const double *src,
struct cfs_tpfa_data *h)
static void
assemble_cell_contrib(struct UnstructuredGrid *G ,
const double *src,
struct cfs_tpfa_data *h )
/* ---------------------------------------------------------------------- */
{
int c1, c2, c, i, f, j1, j2;
int is_neumann;
const double *ctrans = h->pimpl->ctrans;
is_neumann = 1;
for (c = i = 0; c < G->number_of_cells; c++) {
j1 = csrmatrix_elm_index(c, c, h->A);
@ -656,11 +670,6 @@ assemble_cell_contrib(struct UnstructuredGrid *G,
h->A->sa[j1] += ctrans[i];
h->A->sa[j2] -= ctrans[i];
} else if (bc->type[f] == PRESSURE) {
is_neumann = 0;
h->A->sa[j1] += ctrans[i];
h->b [c ] += ctrans[i] * bc->bcval[f];
}
}
@ -669,6 +678,47 @@ assemble_cell_contrib(struct UnstructuredGrid *G,
/* Compressible accumulation term */
h->A->sa[j1] += h->pimpl->accum[c];
}
}
/* ---------------------------------------------------------------------- */
static int
assemble_bc_contrib(struct UnstructuredGrid *G ,
struct FlowBoundaryConditions *fbc,
struct cfs_tpfa_data *h )
/* ---------------------------------------------------------------------- */
{
int c1, c2, c, f, hf;
int is_neumann;
size_t p, j;
const double *ctrans = h->pimpl->ctrans;
is_neumann = 1;
for (p = 0; p < fbc->nbc; p++) {
f = fbc->face[ p ];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
assert ((c1 < 0) ^ (c2 < 0));
c = (c1 >= 0) ? c1 : c2;
if (fbc->type[ p ] == BC_PRESSURE) {
is_neumann = 0;
hf = h->pimpl->f2hf[2*f + (c1 < 0)];
j = csrmatrix_elm_index(c, c, h->A);
h->A->sa[ j ] += ctrans[ hf ];
h->b [ c ] += ctrans[ hf ] * fbc->value[ p ];
}
/* Other boundary condition types not handled. */
}
return is_neumann;
}
@ -735,7 +785,7 @@ assemble_well_contrib(size_t nc,
/* ---------------------------------------------------------------------- */
static void
compute_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *fbc,
int np,
const double *htrans,
const double *pmobf,
@ -746,7 +796,7 @@ compute_fpress(struct UnstructuredGrid *G,
double *scratch_f)
/* ---------------------------------------------------------------------- */
{
int c, i, f, c1, c2;
int c, i, f;
/* Suppress warning about unused parameters. */
(void) np; (void) pmobf; (void) gravcap_f; (void) fflux;
@ -777,12 +827,13 @@ compute_fpress(struct UnstructuredGrid *G,
for (f = 0; f < G->number_of_faces; f++) {
fpress[f] /= scratch_f[f];
}
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == PRESSURE)) {
fpress[f] = bc->bcval[f];
if (fbc != NULL) {
for (i = 0; ((size_t) i) < fbc->nbc; i++) {
if (fbc->type[ i ] == BC_PRESSURE) {
fpress[ fbc->face[ i ] ] = fbc->value[ i ];
}
}
}
}
@ -791,7 +842,7 @@ compute_fpress(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
static void
compute_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *trans,
const double *pmobf,
@ -800,15 +851,14 @@ compute_flux(struct UnstructuredGrid *G,
double *fflux)
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, p;
int f, c1, c2, p, i;
double t, dp, g;
for (f = 0; f < G->number_of_faces; f++) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == FLUX)) {
fflux[f] = bc->bcval[f];
if ((c1 < 0) || (c2 < 0)) {
continue;
}
@ -817,24 +867,42 @@ compute_flux(struct UnstructuredGrid *G,
t += pmobf[f*np + p];
g += pmobf[f*np + p] * gravcap_f[f*np + p];
}
/* t *= trans[f]; */
if ((c1 >= 0) && (c2 >= 0)) {
dp = cpress[c1] - cpress[c2];
} else if (bc->type[f] == PRESSURE) {
if (c1 < 0) {
dp = bc->bcval[f] - cpress[c2];
/* g = -g; */
} else {
dp = cpress[c1] - bc->bcval[f];
}
} else {
/* No BC -> no-flow (== pressure drop offsets gravity) */
dp = -g / t;
}
dp = cpress[c1] - cpress[c2];
fflux[f] = trans[f] * (t*dp + g);
}
if (bc != NULL) {
for (i = 0; ((size_t) i) < bc->nbc; i++) {
f = bc->face[ i ];
if (bc->type[ i ] == BC_FLUX_TOTVOL) {
t = 2*(G->face_cells[2*f + 0] < 0) - 1.0;
fflux[ f ] = t * bc->value[ i ];
}
else if (bc->type[ i ] == BC_PRESSURE) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
t = g = 0.0;
for (p = 0; p < np; p++) {
t += pmobf[f*np + p];
g += pmobf[f*np + p] * gravcap_f[f*np + p];
}
if (c1 < 0) {
dp = bc->value[ i ] - cpress[c1];
}
else {
dp = cpress[c2] - bc->value[ i ];
}
fflux[f] = trans[f] * (t*dp + g);
}
}
}
}
@ -889,6 +957,29 @@ is_incompr(int nc, struct compr_quantities *cq)
}
/* ---------------------------------------------------------------------- */
static void
compute_f2hf_mapping(struct UnstructuredGrid *G ,
struct cfs_tpfa_data *data)
/* ---------------------------------------------------------------------- */
{
int c, nc, i, f, p;
for (i = 0; i < 2 * G->number_of_faces; i++) {
data->pimpl->f2hf[ i ] = -1;
}
for (c = i = 0, nc = G->number_of_cells; c < nc; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[ i ];
p = G->cell_faces[2*f + 0] != c;
data->pimpl->f2hf[2*f + p] = i;
}
}
}
/* ======================================================================
* Public interface below separator.
@ -924,6 +1015,8 @@ cfs_tpfa_construct(struct UnstructuredGrid *G, well_t *W, int nphases)
nwconn = W->well_connpos[ W->number_of_wells ];
}
compute_f2hf_mapping(G, new);
/* Allocate linear system components */
new->b = new->pimpl->ddata + 0;
new->x = new->b + new->A->m;
@ -953,7 +1046,7 @@ void
cfs_tpfa_assemble(struct UnstructuredGrid *G,
double dt,
well_t *W,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
const double *src,
struct compr_quantities *cq,
const double *trans,
@ -973,7 +1066,13 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
compute_psys_contrib(G, W, wdata, bc, cq, dt,
trans, gravcap_f, cpress0, porevol, h);
res_is_neumann = assemble_cell_contrib(G, bc, src, h);
res_is_neumann = 1;
assemble_cell_contrib(G, src, h);
if (bc != NULL) {
res_is_neumann = assemble_bc_contrib(G, bc, h);
}
if ((W != NULL) && (wctrl != NULL)) {
assert (wdata != NULL);
@ -993,7 +1092,7 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_press_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_t *W,
int np,
const double *trans,
@ -1028,7 +1127,7 @@ cfs_tpfa_press_flux(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *htrans,
const double *pmobf,

6
opm/core/pressure/tpfa/cfs_tpfa.h
View File

@ -48,7 +48,7 @@ void
cfs_tpfa_assemble(struct UnstructuredGrid *G,
double dt,
well_t *W,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
const double *src,
struct compr_quantities *cq,
const double *trans,
@ -61,7 +61,7 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
void
cfs_tpfa_press_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_t *W,
int np,
const double *trans,
@ -76,7 +76,7 @@ cfs_tpfa_press_flux(struct UnstructuredGrid *G,
void
cfs_tpfa_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *htrans,
const double *pmobf,