Function cfs_tpfa_residual_assemble() and friends now return singularity flag.

The singularity flag is true if there are no pressure conditions and no
compressibility (so the absolute values of the pressure solution will be
arbitrary).
This commit is contained in:
Atgeirr Flø Rasmussen 2012-08-23 14:00:04 +02:00
parent e96421dbd7
commit 983a29691c
2 changed files with 29 additions and 26 deletions

View File

@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G , cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
double dt , double dt ,
struct cfs_tpfa_res_forces *forces , struct cfs_tpfa_res_forces *forces ,
@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_data *h ) struct cfs_tpfa_res_data *h )
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int res_is_neumann, well_is_neumann, c, np2; int res_is_neumann, well_is_neumann, c, np2, singular;
csrmatrix_zero( h->J); csrmatrix_zero( h->J);
vector_zero (h->J->m, h->F); vector_zero (h->J->m, h->F);
@ -1207,14 +1207,17 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
assemble_sources(dt, forces->src, h); assemble_sources(dt, forces->src, h);
} }
if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) { singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp;
h->J->sa[0] *= 2; if (singular) {
h->J->sa[0] *= 2.0;
} }
return singular;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void int
cfs_tpfa_res_comprock_assemble( cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G , struct UnstructuredGrid *G ,
double dt , double dt ,
@ -1240,29 +1243,18 @@ cfs_tpfa_res_comprock_assemble(
* porevol(pressure)*rock_comp(pressure)/dt * porevol(pressure)*rock_comp(pressure)/dt
*/ */
int c, w, well_is_neumann, rock_is_incomp; int c, rock_is_incomp, singular;
size_t j; size_t j;
double dpv; double dpv;
const struct Wells* W;
/* Assemble usual system (without rock compressibility). */ /* Assemble usual system (without rock compressibility). */
cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
cpress, wpress, porevol0, h); cpress, wpress, porevol0, h);
/* Check if we have only Neumann wells. */
well_is_neumann = 1;
W = forces->wells->W;
for (w = 0; well_is_neumann && w < W->number_of_wells; w++) {
if ((W->ctrls[w]->current >= 0) && /* OPEN? */
(W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) {
well_is_neumann = 0;
}
}
/* If we made a singularity-removing adjustment in the /* If we made a singularity-removing adjustment in the
regular assembly, we undo it here. */ regular assembly, we undo it here. */
if (well_is_neumann && h->pimpl->is_incomp) { if (singular) {
h->J->sa[0] /= 2; h->J->sa[0] /= 2.0;
} }
/* Add new terms to residual and Jacobian. */ /* Add new terms to residual and Jacobian. */
@ -1280,9 +1272,11 @@ cfs_tpfa_res_comprock_assemble(
} }
/* Re-do the singularity-removing adjustment if necessary */ /* Re-do the singularity-removing adjustment if necessary */
if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) { if (rock_is_incomp && singular) {
h->J->sa[0] *= 2; h->J->sa[0] *= 2.0;
} }
return rock_is_incomp && singular;
} }

View File

@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
void void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h); cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
void /* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible and
there are no pressure conditions on wells or boundaries.
Otherwise return 0. */
int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G, cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
double dt, double dt,
struct cfs_tpfa_res_forces *forces, struct cfs_tpfa_res_forces *forces,
@ -72,7 +76,12 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
const double *porevol, const double *porevol,
struct cfs_tpfa_res_data *h); struct cfs_tpfa_res_data *h);
void /* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible, the
rock is incompressible, and there are no pressure conditions on
wells or boundaries.
Otherwise return 0. */
int
cfs_tpfa_res_comprock_assemble( cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G, struct UnstructuredGrid *G,
double dt, double dt,