diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 4f83e42b..317732ee 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -729,9 +729,8 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt, s2 += d2[ p ]; } - /* Negative sign due to completion flux sign convention. */ - h->J->sa[ jc ] -= dt * s1; - h->J->sa[ jw ] -= dt * s2; + h->J->sa[ jc ] += dt * s1; + h->J->sa[ jw ] += dt * s2; } @@ -817,7 +816,7 @@ assemble_completion_to_well(int w, int c, int nc, int np, } -static void +static int assemble_well_contrib(struct cfs_tpfa_res_wells *W , struct compr_quantities_gen *cq , double dt , @@ -826,6 +825,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *W , struct cfs_tpfa_res_data *h ) { int w, i, c, np, np2, nc; + int is_neumann; double pw, dp; double *WI, *gpot, *pmobp; const double *Ac, *dAc; @@ -838,6 +838,8 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *W , gpot = W->data->gpot; pmobp = W->data->phasemob; + is_neumann = 1; + for (w = i = 0; w < W->conn->number_of_wells; w++) { pw = wpress[ w ]; @@ -865,7 +867,12 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *W , if (W->ctrl->ctrl[ w ] != BHP) { h->F[ nc + w ] -= dt * W->ctrl->target[ w ]; } + else { + is_neumann = 0; + } } + + return is_neumann; } @@ -1085,7 +1092,7 @@ cfs_tpfa_res_assemble(grid_t *G, struct cfs_tpfa_res_data *h) /* ---------------------------------------------------------------------- */ { - int res_is_neumann, c, np2; + int res_is_neumann, well_is_neumann, c, np2; csrmatrix_zero( h->J); vector_zero (h->J->m, h->F); @@ -1095,6 +1102,9 @@ cfs_tpfa_res_assemble(grid_t *G, compute_compflux_and_deriv(G, cq->nphases, cpress, trans, cq->phasemobf, gravcap_f, cq->Af, h->pimpl); + res_is_neumann = 1; + well_is_neumann = 1; + np2 = cq->nphases * cq->nphases; for (c = 0; c < G->number_of_cells; c++, zc += cq->nphases) { @@ -1110,7 +1120,8 @@ cfs_tpfa_res_assemble(grid_t *G, compute_well_compflux_and_deriv(forces->W, cq->nphases, cpress, wpress, h->pimpl); - assemble_well_contrib(forces->W, cq, dt, cpress, wpress, h); + well_is_neumann = assemble_well_contrib(forces->W, cq, dt, + cpress, wpress, h); } if ((forces != NULL) && (forces->src != NULL)) { @@ -1118,9 +1129,7 @@ cfs_tpfa_res_assemble(grid_t *G, assemble_sources(dt, forces->src, h); } - res_is_neumann = 1; - - if (res_is_neumann && h->pimpl->is_incomp) { + if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) { h->J->sa[0] *= 2; } } diff --git a/opm/core/utility/cpgpreprocess/cgridinterface.c b/opm/core/utility/cpgpreprocess/cgridinterface.c index 2d5ea66e..772d1eeb 100644 --- a/opm/core/utility/cpgpreprocess/cgridinterface.c +++ b/opm/core/utility/cpgpreprocess/cgridinterface.c @@ -134,8 +134,6 @@ void preprocess (const struct grdecl *in, process_grdecl(in, tol, &pg); - ok = allocate_geometry(base); - /* * General grid interface */ @@ -160,6 +158,7 @@ void preprocess (const struct grdecl *in, base->cell_centroids = NULL; base->cell_volumes = NULL; + ok = allocate_geometry(base); /* * Cornerpoint grid interface