diff --git a/src/TPFACompressiblePressureSolver.hpp b/src/TPFACompressiblePressureSolver.hpp index ea372120..ed5cffb2 100644 --- a/src/TPFACompressiblePressureSolver.hpp +++ b/src/TPFACompressiblePressureSolver.hpp @@ -157,6 +157,7 @@ public: const std::vector& phasemobf, const std::vector& phasemobwellperf, const std::vector& cell_pressure, + const std::vector& gravcapf, const std::vector& wellperf_gpot, const double* surf_dens) { @@ -198,47 +199,15 @@ public: // Assemble the embedded linear system. compr_quantities cq = { 3, &totcompr[0], &voldiscr[0], &cellA[0], &faceA[0], &phasemobf[0] }; - std::vector gravcap_f(3*num_faces, 0.0); - typedef GridAdapter::Vector Vec; - for (int face = 0; face < num_faces; ++face) { - Vec fc = grid_.faceCentroid(face); - for (int local_cell = 0; local_cell < 2; ++local_cell) { - // Total contribution is sum over neighbouring cells. - int cell = grid_.faceCell(face, local_cell); - if (cell == -1) { - // \TODO check that a zero contribution is correct on boundary. - continue; - } - // Compute phase densities in cell. - double phase_dens[3] = { 0.0, 0.0, 0.0 }; - for (int phase = 0; phase < 3; ++phase) { - const double* At = &cellA[9*cell]; // Already transposed since in Fortran order... - for (int comp = 0; comp < 3; ++comp) { - phase_dens[phase] += At[3*phase + comp]*surf_dens[comp]; - } - } - // Compute geometric part. - double gdz = 0.0; - Vec cc = grid_.cellCentroid(cell); - for (int dd = 0; dd < 3; ++dd) { - gdz += (cc[dd] - fc[dd])*gravity_[dd]; - } - if (local_cell == 1) { - gdz *= -1.0; - } - // Add contribution from this cell. - for (int phase = 0; phase < 3; ++phase) { - gravcap_f[3*face + phase] -= gdz*phase_dens[phase]; - } - } - } + + // Call the assembly routine. After this, linearSystem() may be called. cfs_tpfa_assemble(g, dt, wells, &bc, src, - &cq, &trans_[0], &gravcap_f[0], + &cq, &trans_[0], &gravcapf[0], wctrl, wcompl, &cell_pressure[0], &porevol_[0], data_); phasemobf_ = phasemobf; - gravcapf_ = gravcap_f; + gravcapf_ = gravcapf; state_ = Assembled; } @@ -334,7 +303,7 @@ public: data_, &cell_pressures[0], &face_fluxes[0], wpress, wflux); cfs_tpfa_fpress(grid_.c_grid(), &bc, np, &htrans_[0], - &phasemobf_[0], &gravcapf_[0], &cell_pressures[0], + &phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0], &face_fluxes[0], &face_pressures[0]); } diff --git a/src/cfs_tpfa.c b/src/cfs_tpfa.c index 9bba9113..0881159d 100644 --- a/src/cfs_tpfa.c +++ b/src/cfs_tpfa.c @@ -52,6 +52,9 @@ struct cfs_tpfa_impl { double *masstrans_p; /* RB^{-1} [ phase-mobility ] */ double *gravtrans_p; /* RB^{-1} [ grav ] */ + /* Scratch array for face pressure calculation */ + double *scratch_f; + struct densrat_util *ratio; /* Linear storage */ @@ -166,6 +169,8 @@ impl_allocate(grid_t *G, well_t *W, int np) ddata_sz += np * nwperf; /* masstrans_p */ ddata_sz += np * nwperf; /* gravtrans_p */ + ddata_sz += 1 * G->number_of_faces; /* scratch_f */ + new = malloc(1 * sizeof *new); if (new != NULL) { @@ -737,42 +742,45 @@ compute_fpress(grid_t *G, const double *gravcap_f, const double *cpress, const double *fflux, - double *fpress) + double *fpress, + double *scratch_f) /* ---------------------------------------------------------------------- */ { - int c, i, f, p, c1, c2; - double t, s; + int c, i, f, c1, c2; + /* + * Equation used for face pressure pf: + * pf = (h1 p1 + h2 p2) / (h1 + h2) + * where h{12} are half-transmissibilities + * and p{12} are cell pressures. + * + * NOTE: this should be modified to account for gravity and + * Neumann boundaries with nonzero flux! + */ for (f = 0; f < G->number_of_faces; f++) { - fpress[f] = 0.0; + scratch_f[f] = fpress[f] = 0.0; } + /* Temporarily storing (h1 + h2) in scratch[f] + * and (h1 p1 + h2 p2) in fpress[f]. + */ for (c = i = 0; c < G->number_of_cells; c++) { for (; i < G->cell_facepos[c + 1]; i++) { f = G->cell_faces[i]; - - t = 0.0; - for (p = 0; p < np; p++) { - t += pmobf[f*np + p]; - } - t *= htrans[i]; - - s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0; - - fpress[f] += cpress[c] - (s * fflux[f] / t); + scratch_f[f] += htrans[i]; + fpress[f] += htrans[i]*cpress[c]; } } for (f = 0; f < G->number_of_faces; f++) { + fpress[f] /= scratch_f[f]; c1 = G->face_cells[2*f + 0]; c2 = G->face_cells[2*f + 1]; - - fpress[f] /= (c1 >= 0) + (c2 >= 0); - if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == PRESSURE)) { fpress[f] = bc->bcval[f]; } } + } @@ -927,6 +935,9 @@ cfs_tpfa_construct(grid_t *G, well_t *W, int nphases) new->pimpl->wgpot = new->pimpl->wtrans + nwconn; new->pimpl->masstrans_p = new->pimpl->wgpot + nwconn; new->pimpl->gravtrans_p = new->pimpl->masstrans_p + (nphases * nwconn); + + /* Allocate scratch array for face pressure calculations */ + new->pimpl->scratch_f = new->pimpl->gravtrans_p + (nphases * nwconn); } return new; @@ -1012,19 +1023,20 @@ cfs_tpfa_press_flux(grid_t *G, /* ---------------------------------------------------------------------- */ void -cfs_tpfa_fpress(grid_t *G, - flowbc_t *bc, - int np, - const double *htrans, - const double *pmobf, - const double *gravcap_f, - const double *cpress, - const double *fflux, - double *fpress) +cfs_tpfa_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_data *h, + const double *cpress, + const double *fflux, + double *fpress) /* ---------------------------------------------------------------------- */ { compute_fpress(G, bc, np, htrans, pmobf, gravcap_f, - cpress, fflux, fpress); + cpress, fflux, fpress, h->pimpl->scratch_f); } diff --git a/src/cfs_tpfa.h b/src/cfs_tpfa.h index 610870e9..acd809ea 100644 --- a/src/cfs_tpfa.h +++ b/src/cfs_tpfa.h @@ -75,15 +75,16 @@ cfs_tpfa_press_flux(grid_t *G, double *wflux); void -cfs_tpfa_fpress(grid_t *G, - flowbc_t *bc, - int np, - const double *htrans, - const double *pmobf, - const double *gravcap_f, - const double *cpress, - const double *fflux, - double *fpress); +cfs_tpfa_fpress(grid_t *G, + flowbc_t *bc, + int np, + const double *htrans, + const double *pmobf, + const double *gravcap_f, + struct cfs_tpfa_data *h, + const double *cpress, + const double *fflux, + double *fpress); void cfs_tpfa_retrieve_masstrans(grid_t *G,