Enable further compression of boundary conditions.

Specifically, refine the representation to allow a boundary condition
to apply to a set of interfaces rather than just a single face.  This
support is realised by the introduction of new interface
flow_conditions_append_multi().

Update all consumers of struct FlowBoundaryConditions to accommodate
the one-to-many boundary condition interface mapping in the process.
This commit is contained in:
Bård Skaflestad 2012-03-07 19:08:14 +01:00
parent 0d7b7f173f
commit 64f1112997
5 changed files with 163 additions and 86 deletions

View File

@ -19,6 +19,7 @@
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <opm/core/pressure/flow_bc.h>
@ -50,48 +51,63 @@ static void
initialise_structure(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
fbc->nbc = 0;
fbc->cpty = 0;
fbc->nbc = 0;
fbc->cond_cpty = 0;
fbc->face_cpty = 0;
fbc->type = NULL;
fbc->value = NULL;
fbc->face = NULL;
fbc->cond_pos = NULL;
fbc->type = NULL;
fbc->value = NULL;
fbc->face = NULL;
}
/* ---------------------------------------------------------------------- */
static int
expand_tables(size_t nbc,
size_t nf ,
struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
int ok;
int ok_cond, ok_face;
size_t alloc_sz;
void *p1, *p2, *p3;
void *p1, *p2, *p3, *p4;
ok = nbc <= fbc->cpty;
ok_cond = nbc <= fbc->cond_cpty;
ok_face = nf <= fbc->face_cpty;
if (! ok) {
alloc_sz = alloc_size(nbc, fbc->cpty);
if (! ok_cond) {
alloc_sz = alloc_size(nbc, fbc->cond_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 );
p1 = realloc(fbc->type , (alloc_sz + 0) * sizeof *fbc->type );
p2 = realloc(fbc->value , (alloc_sz + 0) * sizeof *fbc->value );
p3 = realloc(fbc->cond_pos, (alloc_sz + 1) * sizeof *fbc->cond_pos);
ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL);
ok_cond = (p1 != NULL) && (p2 != NULL) && (p3 != NULL);
if (ok) {
fbc->type = p1;
fbc->value = p2;
fbc->face = p3;
if (p1 != NULL) { fbc->type = p1; }
if (p2 != NULL) { fbc->value = p2; }
if (p3 != NULL) { fbc->cond_pos = p3; }
fbc->cpty = alloc_sz;
} else {
free(p3); free(p2); free(p1);
if (ok_cond) {
fbc->cond_cpty = alloc_sz;
}
}
return ok;
if (! ok_face) {
alloc_sz = alloc_size(nf, fbc->face_cpty);
p4 = realloc(fbc->face, alloc_sz * sizeof *fbc->face);
ok_face = p4 != NULL;
if (ok_face) {
fbc->face = p4;
fbc->face_cpty = alloc_sz;
}
}
return ok_cond && ok_face;
}
@ -116,12 +132,14 @@ flow_conditions_construct(size_t nbc)
if (fbc != NULL) {
initialise_structure(fbc);
ok = expand_tables(nbc, fbc);
ok = expand_tables(nbc, nbc, fbc);
if (! ok) {
flow_conditions_destroy(fbc);
fbc = NULL;
} else {
fbc->cond_pos[0] = 0;
}
}
@ -138,9 +156,10 @@ flow_conditions_destroy(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
if (fbc != NULL) {
free(fbc->face );
free(fbc->value);
free(fbc->type );
free(fbc->face );
free(fbc->cond_pos);
free(fbc->value );
free(fbc->type );
}
free(fbc);
@ -159,14 +178,38 @@ flow_conditions_append(enum FlowBCType type ,
struct FlowBoundaryConditions *fbc )
/* ---------------------------------------------------------------------- */
{
int ok;
return flow_conditions_append_multi(type, 1, &face, value, fbc);
}
ok = expand_tables(fbc->nbc + 1, fbc);
/* ---------------------------------------------------------------------- */
/* Append a new boundary condition that affects multiple interfaces.
*
* Return one (1) if successful, and zero (0) otherwise. */
/* ---------------------------------------------------------------------- */
int
flow_conditions_append_multi(enum FlowBCType type ,
size_t nfaces,
const int *faces ,
double value ,
struct FlowBoundaryConditions *fbc )
/* ---------------------------------------------------------------------- */
{
int ok;
size_t nbc;
nbc = fbc->nbc;
ok = expand_tables(nbc + 1, fbc->cond_pos[ nbc ] + nfaces, fbc);
if (ok) {
fbc->type [ fbc->nbc ] = type ;
fbc->value[ fbc->nbc ] = value;
fbc->face [ fbc->nbc ] = face ;
memcpy(fbc->face + fbc->cond_pos[ nbc ],
faces, nfaces * sizeof *faces);
fbc->type [ nbc ] = type;
fbc->value[ nbc ] = value;
fbc->cond_pos[ nbc + 1 ] = fbc->cond_pos[ nbc ] + nfaces;
fbc->nbc += 1;
}

View File

@ -38,7 +38,11 @@ enum FlowBCType { BC_NOFLOW ,
* 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 */
size_t face_cpty; /* Internal management. Do not touch */
size_t cond_cpty; /* Internal management. Do not touch */
size_t *cond_pos; /* Indirection pointer into '.face' */
enum FlowBCType *type; /* Condition type */
double *value; /* Condition value (target) */
@ -67,6 +71,16 @@ flow_conditions_append(enum FlowBCType type ,
double value,
struct FlowBoundaryConditions *fbc );
/* Append a new boundary condition that affects multiple interfaces.
*
* Return one (1) if successful, and zero (0) otherwise. */
int
flow_conditions_append_multi(enum FlowBCType type ,
size_t nfaces,
const int *faces ,
double value ,
struct FlowBoundaryConditions *fbc );
/* Clear existing set of boundary conditions */
void

View File

@ -200,17 +200,22 @@ fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int f, i;
int f, i;
size_t j;
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 (fbc->cond_pos[0] == 0);
assert ((0 <= f) && (f < pimpl->nf));
for (i = 0, j = 0; ((size_t) i) < fbc->nbc; i++) {
for (; j < fbc->cond_pos[i + 1]; j++) {
f = fbc->face[ j ];
pimpl->bdry_condition[ f ] = i;
assert ((0 <= f) && (f < pimpl->nf));
pimpl->bdry_condition[ f ] = i;
}
}
}
}

View File

@ -405,7 +405,8 @@ set_dynamic_grav(struct UnstructuredGrid *G ,
struct densrat_util *ratio)
/* ---------------------------------------------------------------------- */
{
int f, p, i, c1, c2;
int f, p, i, c1, c2;
size_t j;
for (f = i = 0; f < G->number_of_faces; f++) {
c1 = G->face_cells[2*f + 0];
@ -425,10 +426,12 @@ 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 (j = bc->cond_pos[i]; j < bc->cond_pos[i + 1]; j++) {
f = bc->face[ j ];
for (p = cq->nphases * f; p < cq->nphases * (f + 1); p++) {
ratio->x[p] = trans[f] * gravcap_f[p] * cq->phasemobf[p];
for (p = cq->nphases * f; p < cq->nphases * (f + 1); p++) {
ratio->x[p] = trans[f] * gravcap_f[p] * cq->phasemobf[p];
}
}
}
}
@ -691,30 +694,32 @@ assemble_bc_contrib(struct UnstructuredGrid *G ,
int c1, c2, c, f, hf;
int is_neumann;
size_t p, j;
size_t p, i, 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;
for (i = fbc->cond_pos[ p ]; i < fbc->cond_pos[p + 1]; i++) {
f = fbc->face[ i ];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
hf = h->pimpl->f2hf[2*f + (c1 < 0)];
assert ((c1 < 0) ^ (c2 < 0));
j = csrmatrix_elm_index(c, c, h->A);
c = (c1 >= 0) ? c1 : c2;
h->A->sa[ j ] += ctrans[ hf ];
h->b [ c ] += ctrans[ hf ] * fbc->value[ p ];
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. */
@ -797,6 +802,7 @@ compute_fpress(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
{
int c, i, f;
size_t j;
/* Suppress warning about unused parameters. */
(void) np; (void) pmobf; (void) gravcap_f; (void) fflux;
@ -832,7 +838,9 @@ compute_fpress(struct UnstructuredGrid *G,
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 ];
for (j = fbc->cond_pos[ i ]; j < fbc->cond_pos[i + 1]; j++) {
fpress[ fbc->face[ j ] ] = fbc->value[ i ];
}
}
}
}
@ -852,6 +860,7 @@ compute_flux(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, p, i;
size_t j;
double t, dp, g;
for (f = 0; f < G->number_of_faces; f++) {
@ -874,32 +883,36 @@ compute_flux(struct UnstructuredGrid *G,
}
if (bc != NULL) {
for (i = 0; ((size_t) i) < bc->nbc; i++) {
f = bc->face[ i ];
for (i = 0, j = 0; ((size_t) i) < bc->nbc; i++) {
for (; j < bc->cond_pos[i + 1]; j++) {
f = bc->face[ j ];
if (bc->type[ i ] == BC_FLUX_TOTVOL) {
t = 2*(G->face_cells[2*f + 0] < 0) - 1.0;
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];
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];
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);
}
if (c1 < 0) {
dp = bc->value[ i ] - cpress[c1];
}
else {
dp = cpress[c2] - bc->value[ i ];
}
fflux[f] = trans[f] * (t*dp + g);
/* Other boundary conditions not handled */
}
}
}

View File

@ -161,7 +161,7 @@ assemble_bc_contrib(struct UnstructuredGrid *G ,
int is_neumann;
int i, f, c1, c2;
size_t j;
size_t j, ix;
is_neumann = 1;
@ -169,17 +169,19 @@ assemble_bc_contrib(struct UnstructuredGrid *G ,
if (bc->type[ i ] == BC_PRESSURE) {
is_neumann = 0;
f = bc->face[ i ];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
for (j = bc->cond_pos[ i ]; j < bc->cond_pos[i + 1]; j++) {
f = bc->face[ j ];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
assert ((c1 < 0) ^ (c2 < 0)); /* BCs on ext. faces only */
assert ((c1 < 0) ^ (c2 < 0)); /* BCs on ext. faces only */
c1 = (c1 >= 0) ? c1 : c2;
j = csrmatrix_elm_index(c1, c1, h->A);
c1 = (c1 >= 0) ? c1 : c2;
ix = csrmatrix_elm_index(c1, c1, h->A);
h->A->sa[ j ] += trans[ f ];
h->b [ c1] += trans[ f ] * bc->value[ i ];
h->A->sa[ ix ] += trans[ f ];
h->b [ c1 ] += trans[ f ] * bc->value[ i ];
}
}
/* Other types currently not handled */