Added gravity parameter, and gravity potential computations.

This commit is contained in:
Atgeirr Flø Rasmussen 2011-04-11 11:04:52 +02:00
parent f70309c76b
commit 4214a6aa1d
2 changed files with 28 additions and 12 deletions

View File

@ -1063,7 +1063,8 @@ cfs_tpfa_impes_maxtime_cell(int c,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_data *h,
const double *dpmobf, const double *dpmobf,
const double *surf_dens) const double *surf_dens,
const double *gravity)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
/* Reference: /* Reference:
@ -1072,10 +1073,11 @@ cfs_tpfa_impes_maxtime_cell(int c,
Capillary pressure parts not included. Capillary pressure parts not included.
*/ */
int i, j, f, c2; int i, j, k, f, c2;
double f11, f12, f21, f22; double f11, f12, f21, f22;
double gsgn, dp, dz, tr, tmob, detF, eqv_flux; double gsgn, dp, dzg, tr, tmob, detF, eqv_flux;
const double *pmob; const double *pmob;
const double *A;
/* This is intended to be compatible with the dune-porsol blackoil /* This is intended to be compatible with the dune-porsol blackoil
code. This library is otherwise not depending on particular code. This library is otherwise not depending on particular
orderings of phases or components, so at some point this orderings of phases or components, so at some point this
@ -1083,14 +1085,12 @@ cfs_tpfa_impes_maxtime_cell(int c,
assert (cq->nphases == 3); assert (cq->nphases == 3);
enum { Water = 0, Oil = 1, Gas = 2 }; enum { Water = 0, Oil = 1, Gas = 2 };
enum { num_phases = 3 }; enum { num_phases = 3 };
double gamma[num_phases]; double rho[num_phases];
double pot[num_phases]; double pot[num_phases];
/* Notation: dpmob[Oil][Water] is d/ds_w(lambda_o) */ /* Notation: dpmob[Oil][Water] is d/ds_w(lambda_o) */
double dpmob[num_phases][num_phases] double dpmob[num_phases][num_phases]
= { {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0} }; = { {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0} };
/* Computing gamma */
/* Filling the dpmob array from available data. /* Filling the dpmob array from available data.
Note that we only need the following combinations: Note that we only need the following combinations:
(Water, Water) (Water, Water)
@ -1120,11 +1120,24 @@ cfs_tpfa_impes_maxtime_cell(int c,
/* Initially only interiour faces */ /* Initially only interiour faces */
if (c2 >= 0) { if (c2 >= 0) {
dp = h->x[c] - h->x[c2]; /* Computing density */
dz = G->cell_centroids[3*c + 2] = G->cell_centroids[3*c2 + 2]; A = cq->Af + f*(cq->nphases)*(cq->nphases);
for (j = 0; j < cq->nphases; ++j) { for (j = 0; j < cq->nphases; ++j) {
pot[j] = fabs(dp - gamma[j]*dz); rho[j] = 0.0;
for (k = 0; k < cq->nphases; ++k) {
rho[j] += A[cq->nphases*j + k]*surf_dens[k];
}
} }
/* Computing gravity potentials */
dp = h->x[c] - h->x[c2];
dzg = 0.0;
for (j = 0; j < G->dimensions; ++j) {
dzg += (G->cell_centroids[G->dimensions*c + j] - G->cell_centroids[G->dimensions*c2 + j])*gravity[j];
}
for (j = 0; j < cq->nphases; ++j) {
pot[j] = fabs(dp - rho[j]*dzg);
}
/* Computing the flux parts f_ij */
pmob = cq->phasemobf + f*cq->nphases; pmob = cq->phasemobf + f*cq->nphases;
tr = trans[f]; tr = trans[f];
tmob = pmob[Water] + pmob[Oil] + pmob[Gas]; tmob = pmob[Water] + pmob[Oil] + pmob[Gas];
@ -1161,14 +1174,16 @@ cfs_tpfa_impes_maxtime(grid_t *G,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_data *h,
const double *dpmobf, const double *dpmobf,
const double *surf_dens) const double *surf_dens,
const double *gravity)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int c; int c;
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, dpmobf, surf_dens); cell_dt = cfs_tpfa_impes_maxtime_cell(c, G, cq, trans, porevol, h,
dpmobf, surf_dens, gravity);
if (cell_dt < max_dt) { if (cell_dt < max_dt) {
max_dt = cell_dt; max_dt = cell_dt;
} }

View File

@ -104,7 +104,8 @@ cfs_tpfa_impes_maxtime(grid_t *G,
const double *porevol, const double *porevol,
struct cfs_tpfa_data *h, struct cfs_tpfa_data *h,
const double *dpmobf, const double *dpmobf,
const double *surf_dens); const double *surf_dens,
const double *gravity);
void void
cfs_tpfa_expl_mass_transport(grid_t *G, cfs_tpfa_expl_mass_transport(grid_t *G,