From 552bfd5f002c69b982056dbe1581807b5f913174 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 9 Nov 2010 12:52:31 +0100 Subject: [PATCH] Add flux and face-pressure computations. Specifically, introduce utilities compute_fpress() and compute_flux(). The former is needed to implement the latter across external boundary faces. Moreover, interface pressure values are needed to evaluate fluid properties on faces (specifically R/B). Add small gateway routine, cfs_tpfa_fpress(), to allow callers to recover interface pressure values. Re-implement cfs_tpfa_press_flux() in terms of compute_fpress() and compute_flux(). Also, add fields 'fpress' and 'fpaccum' to struct cfs_tpfa_impl. --- src/cfs_tpfa.c | 107 ++++++++++++++++++++++++++++++++++++++++++------- src/cfs_tpfa.h | 9 ++++- 2 files changed, 100 insertions(+), 16 deletions(-) diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index f723d0520..1a72b62f8 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -11,7 +11,8 @@ struct cfs_tpfa_impl { - double *fgrav; /* Accumulated grav contrib/face */ + double *fpress; /* Face pressure */ + double *fpaccum; /* Linear storage */ double *ddata; @@ -42,6 +43,7 @@ impl_allocate(grid_t *G) ddata_sz = 2 * G->number_of_cells; /* b, x */ ddata_sz += 1 * G->number_of_faces; /* fgrav */ + ddata_sz += 1 * G->number_of_faces; /* fpaccum */ new = malloc(1 * sizeof *new); @@ -171,6 +173,75 @@ solve_cellsys_core(grid_t *G , } +/* ---------------------------------------------------------------------- */ +static void +compute_fpress(grid_t *G, + flowbc_t *bc, + const double *ctrans, + const double *cpress, + struct cfs_tpfa_data *h) +/* ---------------------------------------------------------------------- */ +{ + int c, i, f; + + for (f = 0; f < G->number_of_faces; f++) { + h->pimpl->fpress [f] = 0.0; + h->pimpl->fpaccum[f] = 0.0; + } + + for (c = i = 0; c < G->number_of_cells; c++) { + for (; i < G->cell_facepos[c + 1]; i++) { + f = G->cell_faces[i]; + + h->pimpl->fpress [f] += ctrans[i] * cpress[c]; + h->pimpl->fpaccum[f] += ctrans[i]; + } + } + + for (f = 0; f < G->number_of_faces; f++) { + h->pimpl->fpress[f] /= h->pimpl->fpaccum[f]; + + if (bc->type[f] == PRESSURE) { + h->pimpl->fpress[f] = bc->bcval[f]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +compute_flux(grid_t *G, + size_t np, + const double *cpress, + const double *ptransf, + struct cfs_tpfa_data *h, + double *fflux) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2, f; + size_t p; + double t, dp; + + for (f = 0; f < G->number_of_faces; f++) { + c1 = G->face_cells[2*f + 0]; + c2 = G->face_cells[2*f + 1]; + + t = 0.0; + for (p = 0; p < np; p++) { t += ptransf[f*np + p]; } + + if ((c1 >= 0) && (c2 >= 0)) { + dp = cpress[c1] - cpress[c2]; + } else if (c1 >= 0) { + dp = cpress[c1] - h->pimpl->fpress[f]; + } else { + dp = h->pimpl->fpress[f] - cpress[c2]; + } + + fflux[f] = t * dp; + } +} + + /* ====================================================================== * Public interface below separator. @@ -197,9 +268,10 @@ cfs_tpfa_construct(grid_t *G) if (new != NULL) { new->b = new->pimpl->ddata; - new->x = new->b + new->A->m; + new->x = new->b + new->A->m; - new->pimpl->fgrav = new->x + new->A->m; + new->pimpl->fpress = new->x + new->A->m; + new->pimpl->fpaccum = new->pimpl->fpress + G->number_of_faces; } return new; @@ -263,27 +335,32 @@ cfs_tpfa_assemble(grid_t *G, /* ---------------------------------------------------------------------- */ void cfs_tpfa_press_flux(grid_t *G, - const double *trans, + size_t np, + flowbc_t *bc, + const double *ctrans, + const double *ptransf, struct cfs_tpfa_data *h, double *cpress, double *fflux) /* ---------------------------------------------------------------------- */ { - int c1, c2, f; - /* Assign cell pressure directly from solution vector */ memcpy(cpress, h->x, G->number_of_cells * sizeof *cpress); - for (f = 0; f < G->number_of_faces; f++) { - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; + compute_fpress(G, bc, ctrans, cpress, h); + compute_flux (G, np, cpress, ptransf, h, fflux); +} - if ((c1 >= 0) && (c2 >= 0)) { - fflux[f] = trans[f] * (cpress[c1] - cpress[c2]); - } else { - fflux[f] = 0.0; - } - } + +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_fpress(grid_t *G, + struct cfs_tpfa_data *h, + double *fpress) +/* ---------------------------------------------------------------------- */ +{ + memcpy(fpress, h->pimpl->fpress, + G->number_of_faces * sizeof *fpress); } diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index 5e3fa63c0..09bfec73f 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -52,10 +52,17 @@ cfs_tpfa_assemble(grid_t *G, void cfs_tpfa_press_flux(grid_t *G, - const double *trans, + size_t np, + flowbc_t *bc, + const double *ctrans, + const double *ptransf, struct cfs_tpfa_data *h, double *cpress, double *fflux); +void +cfs_tpfa_fpress(grid_t *G, + struct cfs_tpfa_data *h, + double *fpress); void cfs_tpfa_destroy(struct cfs_tpfa_data *h);