OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Hook residual formulation up to build.

Browse Source
This commit is contained in:
Bård Skaflestad 2011-10-19 13:01:06 +02:00
parent c8d7ca7fbd
commit 6acd22f7d7
3 changed files with 133 additions and 229 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

266
src/cfs_tpfa_residual.c
View File

@ -8,9 +8,9 @@
#include "flow_bc.h" #include "flow_bc.h"
#include "well.h" #include "well.h"
#include "compr_quant.h" #include "compr_quant_general.h"
#include "trans_tpfa.h" #include "trans_tpfa.h"
#include "cfs_tpfa.h" #include "cfs_tpfa_residual.h"
#include "sparse_sys.h" #include "sparse_sys.h"
@ -35,7 +35,7 @@ struct densrat_util {
}; };
struct cfs_tpfa_impl { struct cfs_tpfa_res_impl {
int is_incomp; int is_incomp;
/* One entry per component per face */ /* One entry per component per face */
@ -110,7 +110,7 @@ allocate_densrat(size_t max_conn, int np)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
static void static void
impl_deallocate(struct cfs_tpfa_impl *pimpl) impl_deallocate(struct cfs_tpfa_res_impl *pimpl)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
if (pimpl != NULL) { if (pimpl != NULL) {
@ -123,12 +123,12 @@ impl_deallocate(struct cfs_tpfa_impl *pimpl)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
static struct cfs_tpfa_impl * static struct cfs_tpfa_res_impl *
impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np) impl_allocate(grid_t *G, size_t max_conn, int np)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
size_t nnu, ngconn, nwperf; size_t nnu, ngconn, nwperf;
struct cfs_tpfa_impl *new; struct cfs_tpfa_res_impl *new;
size_t ddata_sz; size_t ddata_sz;
@ -136,13 +136,8 @@ impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np)
ngconn = G->cell_facepos[ G->number_of_cells ]; ngconn = G->cell_facepos[ G->number_of_cells ];
nwperf = 0; nwperf = 0;
if (W != NULL) {
nnu += W->number_of_wells;
nwperf = W->well_connpos[ W->number_of_wells ];
}
/* Linear system */ /* Linear system */
ddata_sz = 2 * nnu; /* b, x */ ddata_sz = nnu; /* b */
/* Reservoir */ /* Reservoir */
ddata_sz += np * G->number_of_faces ; /* compflux_f */ ddata_sz += np * G->number_of_faces ; /* compflux_f */
@ -170,18 +165,15 @@ impl_allocate(grid_t *G, well_t *W, size_t max_conn, int np)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
static struct CSRMatrix * static struct CSRMatrix *
construct_matrix(grid_t *G, well_t *W) construct_matrix(grid_t *G)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int f, c1, c2, w, i, nc, nnu; int f, c1, c2, i, nc, nnu;
size_t nnz; size_t nnz;
struct CSRMatrix *A; struct CSRMatrix *A;
nc = nnu = G->number_of_cells; nc = nnu = G->number_of_cells;
if (W != NULL) {
nnu += W->number_of_wells;
}
A = csrmatrix_new_count_nnz(nnu); A = csrmatrix_new_count_nnz(nnu);
@ -202,18 +194,6 @@ construct_matrix(grid_t *G, well_t *W)
} }
} }
if (W != NULL) {
/* Well <-> cell connections */
for (w = i = 0; w < W->number_of_wells; w++) {
for (; i < W->well_connpos[w + 1]; i++) {
c1 = W->well_cells[i];
A->ia[ 0 + c1 + 1 ] += 1; /* c -> w */
A->ia[ nc + w + 1 ] += 1; /* w -> c */
}
}
}
nnz = csrmatrix_new_elms_pushback(A); nnz = csrmatrix_new_elms_pushback(A);
if (nnz == 0) { if (nnz == 0) {
csrmatrix_delete(A); csrmatrix_delete(A);
@ -238,18 +218,6 @@ construct_matrix(grid_t *G, well_t *W)
} }
} }
if (W != NULL) {
/* Fill well <-> cell connections */
for (w = i = 0; w < W->number_of_wells; w++) {
for (; i < W->well_connpos[w + 1]; i++) {
c1 = W->well_cells[i];
A->ja[ A->ia[ 0 + c1 + 1 ] ++ ] = nc + w;
A->ja[ A->ia[ nc + w + 1 ] ++ ] = c1 ;
}
}
}
assert ((size_t) A->ia[ nnu ] == nnz); assert ((size_t) A->ia[ nnu ] == nnz);
/* Enforce sorted connection structure per row */ /* Enforce sorted connection structure per row */
@ -358,7 +326,7 @@ compute_compflux_and_deriv(grid_t *G ,
const double *pmobf , const double *pmobf ,
const double *gcapf , const double *gcapf ,
const double *Af , const double *Af ,
struct cfs_tpfa_impl *pimpl ) struct cfs_tpfa_res_impl *pimpl )
{ {
int c1, c2, f, np2; int c1, c2, f, np2;
double dp; double dp;
@ -419,7 +387,7 @@ init_cell_contrib(grid_t *G ,
double pvol , double pvol ,
double dt , double dt ,
const double *z , const double *z ,
struct cfs_tpfa_impl *pimpl) struct cfs_tpfa_res_impl *pimpl)
{ {
int c1, c2, f, i, conn, nconn; int c1, c2, f, i, conn, nconn;
double *cflx, *dcflx; double *cflx, *dcflx;
@ -469,7 +437,7 @@ compute_cell_contrib(grid_t *G ,
const double *z , const double *z ,
const double *Ac , const double *Ac ,
const double *dAc , const double *dAc ,
struct cfs_tpfa_impl *pimpl) struct cfs_tpfa_res_impl *pimpl)
{ {
int c1, c2, f, i, nconn, np2, p; int c1, c2, f, i, nconn, np2, p;
MAT_SIZE_T nrhs; MAT_SIZE_T nrhs;
@ -539,7 +507,7 @@ compute_cell_contrib(grid_t *G ,
static int static int
assemble_cell_contrib(grid_t *G, assemble_cell_contrib(grid_t *G,
int c, int c,
struct cfs_tpfa_data *h) struct cfs_tpfa_res_data *h)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int c1, c2, i, f, j1, j2; int c1, c2, i, f, j1, j2;
@ -547,9 +515,9 @@ assemble_cell_contrib(grid_t *G,
is_neumann = 1; is_neumann = 1;
j1 = csrmatrix_elm_index(c, c, h->A); j1 = csrmatrix_elm_index(c, c, h->J);
h->A->sa[j1] += h->pimpl->ratio->mat_row[ 0 ]; h->J->sa[j1] += h->pimpl->ratio->mat_row[ 0 ];
for (; i < G->cell_facepos[c + 1]; i++) { for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i]; f = G->cell_faces[i];
@ -560,13 +528,13 @@ assemble_cell_contrib(grid_t *G,
c2 = (c1 == c) ? c2 : c1; c2 = (c1 == c) ? c2 : c1;
if (c2 >= 0) { if (c2 >= 0) {
j2 = csrmatrix_elm_index(c, c2, h->A); j2 = csrmatrix_elm_index(c, c2, h->J);
h->A->sa[j2] -= h->pimpl->ratio->mat_row[ i ]; h->J->sa[j2] -= h->pimpl->ratio->mat_row[ i ];
} }
} }
h->b[ c ] = h->pimpl->ratio->residual; h->F[ c ] = h->pimpl->ratio->residual;
return 0; return 0;
} }
@ -713,74 +681,83 @@ is_incompr(int nc, const double *accum)
* Public interface below separator. * Public interface below separator.
* ====================================================================== */ * ====================================================================== */
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
struct cfs_tpfa_data * void
cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
size_t nc, nf, ngconn, nwconn; if (h != NULL) {
struct cfs_tpfa_data *new; csrmatrix_delete(h->J);
impl_deallocate (h->pimpl);
}
new = malloc(1 * sizeof *new); free(h);
}
if (new != NULL) {
new->pimpl = impl_allocate(G, W, maxconn(G), nphases);
new->A = construct_matrix(G, W);
if ((new->pimpl == NULL) || (new->A == NULL)) { /* ---------------------------------------------------------------------- */
cfs_tpfa_destroy(new); struct cfs_tpfa_res_data *
new = NULL; cfs_tpfa_res_construct(grid_t *G, int nphases)
/* ---------------------------------------------------------------------- */
{
size_t nc, nf, ngconn;
struct cfs_tpfa_res_data *h;
h = malloc(1 * sizeof *h);
if (h != NULL) {
h->pimpl = impl_allocate(G, maxconn(G), nphases);
h->J = construct_matrix(G);
if ((h->pimpl == NULL) || (h->J == NULL)) {
cfs_tpfa_res_destroy(h);
h = NULL;
} }
} }
if (new != NULL) { if (h != NULL) {
nc = G->number_of_cells; nc = G->number_of_cells;
nf = G->number_of_faces; nf = G->number_of_faces;
ngconn = G->cell_facepos[nc]; ngconn = G->cell_facepos[nc];
nwconn = 0;
if (W != NULL) {
nwconn = W->well_connpos[ W->number_of_wells ];
}
/* Allocate linear system components */ /* Allocate linear system components */
new->b = new->pimpl->ddata + 0; h->F = h->pimpl->ddata + 0;
new->x = new->b + new->A->m;
new->pimpl->compflux_f = new->x + new->A->m; h->pimpl->compflux_f = h->F + h->J->m;
new->pimpl->compflux_deriv_f = h->pimpl->compflux_deriv_f =
new->pimpl->compflux_f + (nphases * nf); h->pimpl->compflux_f + (nphases * nf);
new->pimpl->flux_work = h->pimpl->flux_work =
new->pimpl->compflux_deriv_f + (nphases * 2 * nf); h->pimpl->compflux_deriv_f + (nphases * 2 * nf);
new->pimpl->scratch_f = h->pimpl->scratch_f =
new->pimpl->flux_work + (nphases * (1 + 2)); h->pimpl->flux_work + (nphases * (1 + 2));
} }
return new; return h;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void void
cfs_tpfa_assemble(grid_t *G, cfs_tpfa_res_assemble(grid_t *G,
double dt, double dt,
flowbc_t *bc, flowbc_t *bc,
const double *src, const double *src,
const double *zc, const double *zc,
struct compr_quantities *cq, struct compr_quantities_gen *cq,
const double *trans, const double *trans,
const double *gravcap_f, const double *gravcap_f,
const double *cpress, const double *cpress,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h) struct cfs_tpfa_res_data *h)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int res_is_neumann, c, np2; int res_is_neumann, c, np2;
csrmatrix_zero( h->A); csrmatrix_zero( h->J);
vector_zero (h->A->m, h->b); vector_zero (h->J->m, h->F);
h->pimpl->is_incomp = 1; h->pimpl->is_incomp = 1;
@ -814,70 +791,37 @@ cfs_tpfa_assemble(grid_t *G,
if (res_is_neumann && h->pimpl->is_incomp) { if (res_is_neumann && h->pimpl->is_incomp) {
/*&& well_is_neumann && h->pimpl->is_incomp) {*/ /*&& well_is_neumann && h->pimpl->is_incomp) {*/
h->A->sa[0] *= 2; h->J->sa[0] *= 2;
} }
} }
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_res_flux(grid_t *G,
flowbc_t *bc,
int np,
const double *trans,
const double *pmobf,
const double *gravcap_f,
const double *cpress,
double *fflux)
/* ---------------------------------------------------------------------- */
{
compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux);
}
#if 0 #if 0
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void void
cfs_tpfa_press_increment(grid_t *G, cfs_tpfa_res_fpress(grid_t *G,
well_t *W,
struct cfs_tpfa_data *h,
double *cpress_inc,
double *wpress_inc)
/* ---------------------------------------------------------------------- */
{
/* Assign cell pressure directly from solution vector */
memcpy(cpress_inc, h->x, G->number_of_cells * sizeof *cpress_inc);
if (W != NULL) {
assert (wpress_inc != NULL);
/* Assign well BHP directly from solution vector */
memcpy(wpress_inc, h->x + G->number_of_cells,
W->number_of_wells * sizeof *wpress_inc);
}
}
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_flux(grid_t *G,
flowbc_t *bc,
well_t *W,
int np,
const double *trans,
const double *pmobf,
const double *gravcap_f,
const double *cpress,
const double *wpress,
struct completion_data *wdata,
double *fflux,
double *wflux)
/* ---------------------------------------------------------------------- */
{
compute_flux(G, bc, np, trans, pmobf, gravcap_f, cpress, fflux);
if ((W != NULL) && (wdata != NULL)) {
assert (wpress != NULL);
assert (wflux != NULL);
compute_wflux(W, np, wdata, cpress, wpress, wflux);
}
}
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_fpress(grid_t *G,
flowbc_t *bc, flowbc_t *bc,
int np, int np,
const double *htrans, const double *htrans,
const double *pmobf, const double *pmobf,
const double *gravcap_f, const double *gravcap_f,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
const double *cpress, const double *cpress,
const double *fflux, const double *fflux,
double *fpress) double *fpress)
@ -890,9 +834,9 @@ cfs_tpfa_fpress(grid_t *G,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void void
cfs_tpfa_retrieve_masstrans(grid_t *G, cfs_tpfa_res_retrieve_masstrans(grid_t *G,
int np, int np,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
double *masstrans_f) double *masstrans_f)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
@ -903,9 +847,9 @@ cfs_tpfa_retrieve_masstrans(grid_t *G,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void void
cfs_tpfa_retrieve_gravtrans(grid_t *G, cfs_tpfa_res_retrieve_gravtrans(grid_t *G,
int np, int np,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
double *gravtrans_f) double *gravtrans_f)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
@ -916,12 +860,12 @@ cfs_tpfa_retrieve_gravtrans(grid_t *G,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
static double static double
cfs_tpfa_impes_maxtime_cell(int c, cfs_tpfa_res_impes_maxtime_cell(int c,
grid_t *G, grid_t *G,
struct compr_quantities *cq, struct compr_quantities *cq,
const double *trans, const double *trans,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
const double *dpmobf, const double *dpmobf,
const double *surf_dens, const double *surf_dens,
const double *gravity) const double *gravity)
@ -1033,11 +977,11 @@ cfs_tpfa_impes_maxtime_cell(int c,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double double
cfs_tpfa_impes_maxtime(grid_t *G, cfs_tpfa_res_impes_maxtime(grid_t *G,
struct compr_quantities *cq, struct compr_quantities *cq,
const double *trans, const double *trans,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
const double *dpmobf, const double *dpmobf,
const double *surf_dens, const double *surf_dens,
const double *gravity) const double *gravity)
@ -1047,7 +991,7 @@ cfs_tpfa_impes_maxtime(grid_t *G,
double max_dt, cell_dt; double max_dt, cell_dt;
max_dt = 1e100; max_dt = 1e100;
for (c = 0; c < G->number_of_cells; ++c) { for (c = 0; c < G->number_of_cells; ++c) {
cell_dt = cfs_tpfa_impes_maxtime_cell(c, G, cq, trans, porevol, h, cell_dt = cfs_tpfa_res_impes_maxtime_cell(c, G, cq, trans, porevol, h,
dpmobf, surf_dens, gravity); dpmobf, surf_dens, gravity);
if (cell_dt < max_dt) { if (cell_dt < max_dt) {
max_dt = cell_dt; max_dt = cell_dt;
@ -1060,13 +1004,13 @@ cfs_tpfa_impes_maxtime(grid_t *G,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void void
cfs_tpfa_expl_mass_transport(grid_t *G, cfs_tpfa_res_expl_mass_transport(grid_t *G,
well_t *W, well_t *W,
struct completion_data *wdata, struct completion_data *wdata,
int np, int np,
double dt, double dt,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_res_data *h,
double *surf_vol) double *surf_vol)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
@ -1130,17 +1074,3 @@ cfs_tpfa_expl_mass_transport(grid_t *G,
} }
} }
#endif #endif
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_destroy(struct cfs_tpfa_data *h)
/* ---------------------------------------------------------------------- */
{
if (h != NULL) {
csrmatrix_delete(h->A);
impl_deallocate (h->pimpl);
}
free(h);
}

67
src/cfs_tpfa_residual.h
View File

@ -28,52 +28,48 @@
extern "C" { extern "C" {
#endif #endif
struct cfs_tpfa_impl; struct cfs_tpfa_res_impl;
struct CSRMatrix; struct CSRMatrix;
struct compr_quantities; struct compr_quantities_gen;
struct cfs_tpfa_data { struct cfs_tpfa_res_data {
struct CSRMatrix *A; struct CSRMatrix *J;
double *b; double *F;
double *x;
struct cfs_tpfa_impl *pimpl; struct cfs_tpfa_res_impl *pimpl;
}; };
struct cfs_tpfa_data * struct cfs_tpfa_res_data *
cfs_tpfa_construct(grid_t *G, well_t *W, int nphases); cfs_tpfa_res_construct(grid_t *G, int nphases);
void void
cfs_tpfa_assemble(grid_t *G, cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
double dt,
flowbc_t *bc,
const double *src,
const double *zc,
struct compr_quantities *cq,
const double *trans,
const double *gravcap_f,
const double *cpress,
const double *porevol,
struct cfs_tpfa_data *h);
void void
cfs_tpfa_press_increment(grid_t *G, cfs_tpfa_res_assemble(grid_t *G,
struct cfs_tpfa_data *h, double dt,
double *cpress_inc); flowbc_t *bc,
const double *src,
const double *zc,
struct compr_quantities_gen *cq,
const double *trans,
const double *gravcap_f,
const double *cpress,
const double *porevol,
struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_flux(grid_t *G,
flowbc_t *bc,
int np,
const double *trans,
const double *pmobf,
const double *gravcap_f,
const double *cpress,
double *fflux);
#if 0 #if 0
void
cfs_tpfa_flux(grid_t *G,
flowbc_t *bc,
int np,
const double *trans,
const double *pmobf,
const double *gravcap_f,
const double *cpress,
double *fflux);
#endif
void void
cfs_tpfa_fpress(grid_t *G, cfs_tpfa_fpress(grid_t *G,
flowbc_t *bc, flowbc_t *bc,
@ -118,8 +114,7 @@ cfs_tpfa_expl_mass_transport(grid_t *G,
struct cfs_tpfa_data *h, struct cfs_tpfa_data *h,
double *surf_vol); double *surf_vol);
void #endif
cfs_tpfa_destroy(struct cfs_tpfa_data *h);
#ifdef __cplusplus #ifdef __cplusplus
} }

29
src/compr_quant_general.h
View File

@ -22,13 +22,11 @@
#include <stddef.h> #include <stddef.h>
#include "grid.h"
#ifdef __cplusplus #ifdef __cplusplus
extern "C" { extern "C" {
#endif #endif
struct compr_quantities { struct compr_quantities_gen {
int nphases; /* Number of phases/components */ int nphases; /* Number of phases/components */
double *Ac; /* RB^{-1} per cell */ double *Ac; /* RB^{-1} per cell */
@ -38,30 +36,11 @@ struct compr_quantities {
double *voldiscr; /* Volume discrepancy per cell */ double *voldiscr; /* Volume discrepancy per cell */
}; };
struct compr_quantities * struct compr_quantities_gen *
compr_quantities_allocate(size_t nc, size_t nf, int np); compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
void void
compr_quantities_deallocate(struct compr_quantities *cq); compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);
void
compr_flux_term(grid_t *G,
const double *fflux,
const double *zeta,
double *Biv);
void
compr_accum_term(size_t nc,
double dt,
const double *porevol,
const double *totcompr,
double *P);
void
compr_src_add_press_accum(size_t nc,
const double *p0,
const double *P,
double *src);
#ifdef __cplusplus #ifdef __cplusplus
} }