From c723f7d4f75e40f884f90565c8e36009fdb269fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 26 Aug 2011 15:30:20 +0200 Subject: [PATCH] Add a simple QFS for testing from C. Not integrated into Autotools build system. Compile as gcc -g -Wall -ansi -pedantic -Wextra test_cfs_tpfa.c cfs_tpfa.c \ well.c flow_bc.c trans_tpfa.c sparse_sys.c compr_quant.c \ -lumfpack -llapack -lblas -lm (or variants thereof). --- src/compr_quant.h | 10 +- src/test_cfs_tpfa.c | 690 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 695 insertions(+), 5 deletions(-) create mode 100644 src/test_cfs_tpfa.c diff --git a/src/compr_quant.h b/src/compr_quant.h index 01cf2f4b..cacfb551 100644 --- a/src/compr_quant.h +++ b/src/compr_quant.h @@ -31,11 +31,11 @@ extern "C" { struct compr_quantities { int nphases; /* Number of phases/components */ - const double *totcompr; /* Total compressibility per cell */ - const double *voldiscr; /* Volume discrepancy per cell */ - const double *Ac; /* RB^{-1} per cell */ - const double *Af; /* RB^{-1} per face */ - const double *phasemobf; /* Phase mobility per face */ + double *totcompr; /* Total compressibility per cell */ + double *voldiscr; /* Volume discrepancy per cell */ + double *Ac; /* RB^{-1} per cell */ + double *Af; /* RB^{-1} per face */ + double *phasemobf; /* Phase mobility per face */ }; void diff --git a/src/test_cfs_tpfa.c b/src/test_cfs_tpfa.c new file mode 100644 index 00000000..6becf0ae --- /dev/null +++ b/src/test_cfs_tpfa.c @@ -0,0 +1,690 @@ +#include +#include +#include +#include + +#include + +#include "cfs_tpfa.h" +#include "sparse_sys.h" +#include "compr_quant.h" +#include "trans_tpfa.h" + +#include "grid.h" +#include "flow_bc.h" +#include "well.h" + + +struct CSCMatrix { + UF_long n; + UF_long nnz; + + UF_long *p; + UF_long *i; + double *x; +}; + + +/* ---------------------------------------------------------------------- */ +static void +csc_deallocate(struct CSCMatrix *csc) +/* ---------------------------------------------------------------------- */ +{ + if (csc != NULL) { + if (csc->x != NULL) { free(csc->x); } + if (csc->i != NULL) { free(csc->i); } + if (csc->p != NULL) { free(csc->p); } + + free(csc); + } +} + + +/* ---------------------------------------------------------------------- */ +static struct CSCMatrix * +csc_allocate(UF_long n, UF_long nnz) +/* ---------------------------------------------------------------------- */ +{ + struct CSCMatrix *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->p = malloc((n + 1) * sizeof *new->p); + new->i = malloc(nnz * sizeof *new->i); + new->x = malloc(nnz * sizeof *new->x); + + if ((new->p == NULL) || (new->i == NULL) || (new->x == NULL)) { + csc_deallocate(new); + new = NULL; + } else { + new->n = n; + new->nnz = nnz; + } + } + + return new; +} + + +/* ---------------------------------------------------------------------- */ +static void +csr_to_csc(const int *ia, + const int *ja, + const double *sa, + struct CSCMatrix *csc) +/* ---------------------------------------------------------------------- */ +{ + UF_long i, nz; + + /* Clear garbage, prepare for counting */ + for (i = 0; i <= csc->n; i++) { csc->p[i] = 0; } + + /* Count column connections */ + for (nz = 0; nz < csc->nnz; nz++) { + csc->p[ ja[nz] + 1 ] += 1; + } + + /* Define column start pointers */ + for (i = 1; i <= csc->n; i++) { + csc->p[0] += csc->p[i]; + csc->p[i] = csc->p[0] - csc->p[i]; + } + + assert (csc->p[0] == csc->nnz); + + /* Fill matrix whilst defining column end pointers */ + for (i = nz = 0; i < csc->n; i++) { + for (; nz < ia[i + 1]; nz++) { + csc->i[ csc->p[ ja[nz] + 1 ] ] = i; /* Insertion sort */ + csc->x[ csc->p[ ja[nz] + 1 ] ] = sa[nz]; /* Insert mat elem */ + + csc->p [ ja[nz] + 1 ] += 1; /* Advance col ptr */ + } + } + + assert (csc->p[csc->n] == csc->nnz); + + csc->p[0] = 0; +} + + +/* ---------------------------------------------------------------------- */ +static void +solve_umfpack(struct CSCMatrix *csc, const double *b, double *x) +/* ---------------------------------------------------------------------- */ +{ + void *Symbolic, *Numeric; + double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL]; + + umfpack_dl_defaults(Control); + + umfpack_dl_symbolic(csc->n, csc->n, csc->p, csc->i, csc->x, + &Symbolic, Control, Info); + umfpack_dl_numeric (csc->p, csc->i, csc->x, + Symbolic, &Numeric, Control, Info); + + umfpack_dl_free_symbolic(&Symbolic); + + umfpack_dl_solve(UMFPACK_A, csc->p, csc->i, csc->x, x, b, + Numeric, Control, Info); + + umfpack_dl_free_numeric(&Numeric); +} + + +/*---------------------------------------------------------------------------*/ +void +call_UMFPACK(struct CSRMatrix *A, double *b, double *x) +/*---------------------------------------------------------------------------*/ +{ + struct CSCMatrix *csc; + + csc = csc_allocate(A->m, A->ia[A->m]); + + if (csc != NULL) { + csr_to_csc(A->ia, A->ja, A->sa, csc); + + solve_umfpack(csc, b, x); + } + + csc_deallocate(csc); +} + +static void +deallocate_cart_grid(grid_t *G) +{ + if (G != NULL) { + free(G->node_coordinates); + + free(G->face_nodes); + free(G->face_nodepos); + free(G->face_cells); + free(G->face_centroids); + free(G->face_normals); + free(G->face_areas); + + free(G->cell_faces); + free(G->cell_facepos); + free(G->cell_centroids); + free(G->cell_volumes); + } + + free(G); +} + +grid_t * +cart_grid(int nx, int ny, int nz) +{ + /* The following applies to any grid derived from base_grid_t. */ + int i,j,k; + int Nx = nx+1; + int Ny = ny+1; + int Nz = nz+1; + + grid_t *G = malloc(1 * sizeof *G); + + G->dimensions = 3; + + int nxf = Nx*ny*nz; + int nyf = nx*Ny*nz; + int nzf = nx*ny*Nz; + + G->number_of_cells = nx*ny*nz; + G->number_of_faces = nxf+nyf+nzf; + G->number_of_nodes = Nx*Ny*Nz; + + G->node_coordinates = malloc(G->number_of_nodes * 3 * sizeof *(G->node_coordinates)); + + G->face_nodes = malloc(G->number_of_faces * 4 * sizeof *(G->face_nodes)); + G->face_nodepos = malloc((G->number_of_faces+1) * sizeof *(G->face_nodepos)); + G->face_cells = malloc(G->number_of_faces * 2 * sizeof *(G->face_cells)); + G->face_centroids = malloc(G->number_of_faces * 3 * sizeof *(G->face_centroids)); + G->face_normals = malloc(G->number_of_faces * 3 * sizeof *(G->face_normals)); + G->face_areas = malloc(G->number_of_faces * 1 * sizeof *(G->face_areas)); + + G->cell_faces = malloc(G->number_of_cells * 6 * sizeof *(G->cell_faces)); + G->cell_facepos = malloc((G->number_of_cells+1) * sizeof *(G->cell_facepos)); + G->cell_centroids = malloc(G->number_of_cells * 3 * sizeof *(G->cell_centroids)); + G->cell_volumes = malloc(G->number_of_cells * 1 * sizeof *(G->cell_volumes)); + + + int *cfaces = G->cell_faces; + int *cfacepos = G->cell_facepos; + double *ccentroids = G->cell_centroids; + double *cvolumes = G->cell_volumes; + for (k=0; kface_nodes; + int *fnodepos = G->face_nodepos; + int *fcells = G->face_cells; + double *fnormals = G->face_normals; + double *fcentroids = G->face_centroids; + double *fareas = G->face_areas; + + /* Faces with x-normal */ + for (k=0; knode_coordinates; + for (k=0; kphasemobf); + free(cq->Af); + free(cq->Ac); + free(cq->voldiscr); + free(cq->totcompr); + } + + free(cq); +} + +static struct compr_quantities * +allocate_cq(size_t nc, size_t nf, int np) +{ + struct compr_quantities *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->totcompr = malloc(nc * sizeof *new->totcompr ); + new->voldiscr = malloc(nc * sizeof *new->voldiscr ); + new->Ac = malloc(nc * np * np * sizeof *new->Ac ); + new->Af = malloc(nf * np * np * sizeof *new->Af ); + new->phasemobf = malloc(nf * np * sizeof *new->phasemobf); + + if ((new->totcompr == NULL) || (new->voldiscr == NULL) || + (new->Ac == NULL) || (new->Af == NULL) || + (new->phasemobf == NULL)) { + deallocate_cq(new); + new = NULL; + } else { + new->nphases = np; + } + } + + return new; +} + +static void +vector_assign(size_t n, double value, double *v) +{ + size_t i; + + for (i = 0; i < n; i++) { + v[i] = value; + } +} + +static void +vector_ones(size_t n, double *v) +{ + vector_assign(n, 1.0, v); +} + +static void +set_incompressible(size_t nc, size_t nf, struct compr_quantities *cq) +{ + size_t i, j, np; + + np = cq->nphases; + + vector_zero(nc , cq->totcompr ); + vector_zero(nc , cq->voldiscr ); + vector_zero(nc * np * np, cq->Ac); + vector_zero(nf * np * np, cq->Af); + vector_ones(nf * np , cq->phasemobf); + + for (i = 0; i < nc; i++) { + for (j = 0; j < np; j++) { + cq->Ac[i*np*np + j*(np + 1)] = 1.0; + } + } + + for (i = 0; i < nf; i++) { + for (j = 0; j < np; j++) { + cq->Af[i*np*np + j*(np + 1)] = 1.0; + } + } +} + + +static void +set_homoperm(size_t nc, size_t dim, double *perm) +{ + size_t i, j; + + vector_zero(nc * dim * dim, perm); + + for (i = 0; i < nc; i++) { + for (j = 0; j < dim; j++) { + perm[i*dim*dim + j*(dim + 1)] = 1.0; + } + } +} + + +static void +destroy_wells(well_t *W) +{ + if (W != NULL) { + free(W->well_cells); + free(W->well_connpos); + } + + free(W); +} + + +static well_t * +create_wells(grid_t *G) +{ + well_t *new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->well_connpos = malloc((2 + 1) * sizeof *new->well_connpos); + new->well_cells = malloc(2 * sizeof *new->well_cells ); + + if ((new->well_connpos == NULL) || + (new->well_cells == NULL)) { + destroy_wells(new); + new = NULL; + } else { + new->number_of_wells = 2; + + new->well_connpos[0] = 0; + new->well_connpos[1] = 1; + new->well_connpos[2] = 2; + + new->well_cells[0] = 0; + new->well_cells[1] = G->number_of_cells - 1; + } + } + + return new; +} + + +static void +destroy_well_control(well_control_t *ctrl) +{ + if (ctrl != NULL) { + free(ctrl->target); + free(ctrl->ctrl); + free(ctrl->type); + } + + free(ctrl); +} + + +static well_control_t * +create_well_control(void) +{ + well_control_t *new; + + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->type = malloc(2 * sizeof *new->type ); + new->ctrl = malloc(2 * sizeof *new->ctrl ); + new->target = malloc(2 * sizeof *new->target); + + if ((new->type == NULL) || (new->ctrl == NULL) || + (new->target == NULL)) { + destroy_well_control(new); + new = NULL; + } else { + new->type[0] = INJECTOR; new->type[1] = PRODUCER; + new->ctrl[0] = new->ctrl[1] = BHP; + new->target[0] = 2; + new->target[1] = 1; + } + } + + return new; +} + + +static void +destroy_completion_data(struct completion_data *cd) +{ + if (cd != NULL) { + free(cd->WI); + free(cd->gpot); + free(cd->A); + free(cd->phasemob); + } + + free(cd); +} + + +static struct completion_data * +create_completion_data(well_t *W, int np) +{ + size_t i, j, totconn; + struct completion_data *new; + + totconn = W->well_connpos[ W->number_of_wells ]; + new = malloc(1 * sizeof *new); + + if (new != NULL) { + new->WI = malloc(totconn * sizeof *new->WI); + new->gpot = malloc(totconn * np * sizeof *new->gpot); + new->A = malloc(totconn * np * np * sizeof *new->A); + new->phasemob = malloc(totconn * np * sizeof *new->phasemob); + + if ((new->WI == NULL) || + (new->gpot == NULL) || + (new->A == NULL) || + (new->phasemob == NULL)) { + destroy_completion_data(new); + new = NULL; + } else { + vector_ones(totconn , new->WI); + vector_zero(totconn * np , new->gpot); + vector_zero(totconn * np * np, new->A); + vector_ones(totconn * np , new->phasemob); + + for (i = 0; i < totconn; i++) { + for (j = 0; j < (size_t)np; j++) { + new->A[i*np*np + j*(np + 1)] = 1.0; + } + } + } + } + + return new; +} + + +int +main(void) +{ + int i, nphases, dim; + flowbc_t *bc; + double dt; + double *htrans, *trans, *perm, *totmob, *src; + double *gravcap_f, *porevol, *cpress0; + double *cpress, *fpress, *fflux, *wpress, *wflux; + + grid_t *G; + well_t *W; + well_control_t *wctrl; + struct completion_data *wdata; + + struct compr_quantities *cq; + + struct cfs_tpfa_data *h; + + nphases = 2; + + G = cart_grid(5, 5, 1); + W = create_wells(G); + wctrl = create_well_control(); + wdata = create_completion_data(W, nphases); + + dim = G->dimensions; + + htrans = malloc(G->cell_facepos[ G->number_of_cells ] * sizeof *htrans); + trans = malloc(G->number_of_faces * sizeof *trans ); + perm = malloc(dim * dim * G->number_of_cells * sizeof *perm ); + totmob = malloc(G->number_of_cells * sizeof *totmob); + src = malloc(G->number_of_cells * sizeof *src ); + + gravcap_f = malloc(G->number_of_faces * nphases * sizeof *gravcap_f); + porevol = malloc(G->number_of_cells * sizeof *porevol ); + cpress0 = malloc(G->number_of_cells * sizeof *cpress0 ); + cpress = malloc(G->number_of_cells * sizeof *cpress ); + fpress = malloc(G->number_of_faces * sizeof *fpress ); + fflux = malloc(G->number_of_faces * sizeof *fflux ); + + wpress = malloc(W->number_of_wells * sizeof *wpress); + wflux = malloc(W->well_connpos[ W->number_of_wells ] * sizeof *wflux); + + bc = allocate_flowbc(G->number_of_faces); + + set_homoperm(G->number_of_cells, G->dimensions, perm); + + cq = allocate_cq(G->number_of_cells, G->number_of_faces, nphases); + set_incompressible(G->number_of_cells, G->number_of_faces, cq); + + h = cfs_tpfa_construct(G, W, nphases); + + vector_assign(G->number_of_cells, nphases, totmob); + + tpfa_htrans_compute(G, perm, htrans); + tpfa_eff_trans_compute(G, totmob, htrans, trans); + + vector_zero(G->number_of_cells, src); + + dt = 1; + cfs_tpfa_assemble(G, dt, W, bc, src, cq, trans, gravcap_f, wctrl, wdata, + cpress0, porevol, + h); + + call_UMFPACK(h->A, h->b, h->x); + + cfs_tpfa_press_flux(G, bc, W, nphases, trans, cq->phasemobf, gravcap_f, + wdata, h, cpress, fflux, wpress, wflux); + + cfs_tpfa_fpress(G, bc, nphases, htrans, cq->phasemobf, gravcap_f, h, + cpress, fflux, fpress); + + for (i = 0; i < G->number_of_cells; i++) { + fprintf(stderr, "press(%02d) = %g;\n", i + 1, cpress[i]); + } + + cfs_tpfa_destroy(h); + deallocate_cq(cq); + deallocate_flowbc(bc); + + free(wflux); free(wpress); + free(fflux); free(fpress); free(cpress); + free(cpress0); free(porevol); free(gravcap_f); + free(src); free(totmob); free(perm); free(trans); free(htrans); + + deallocate_cart_grid(G); + + return 0; +}