diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 6568e8cc3..02f675567 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -87,7 +87,7 @@ namespace Opm double perm_threshold) { const int dim = 3; - const int num_global_cells = numGlobalCells(parser); + const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; const int nc = grid.number_of_cells; ASSERT (num_global_cells > 0); diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 154f06336..d28db5a33 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -58,6 +59,7 @@ namespace Opm /// to change. CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -66,6 +68,7 @@ namespace Opm const struct Wells* wells) : grid_(grid), props_(props), + rock_comp_props_(rock_comp_props), linsolver_(linsolver), residual_tol_(residual_tol), change_tol_(change_tol), @@ -74,7 +77,6 @@ namespace Opm wells_(wells), htrans_(grid.cell_facepos[ grid.number_of_cells ]), trans_ (grid.number_of_faces), - porevol_(grid.number_of_cells), allcells_(grid.number_of_cells) { if (wells_ && (wells_->number_of_phases != props.numPhases())) { @@ -86,7 +88,12 @@ namespace Opm UnstructuredGrid* gg = const_cast(&grid_); tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); tpfa_trans_compute(gg, &htrans_[0], &trans_[0]); - computePorevolume(grid_, props.porosity(), porevol_); + // If we have rock compressibility, pore volumes are updated + // in the compute*() methods, otherwise they are constant and + // hence may be computed here. + if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { + computePorevolume(grid_, props.porosity(), porevol_); + } for (int c = 0; c < grid.number_of_cells; ++c) { allcells_[c] = c; } @@ -230,6 +237,9 @@ namespace Opm const WellState& /*well_state*/) { computeWellPotentials(state); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_); + } } @@ -252,6 +262,8 @@ namespace Opm // std::vector face_gravcap_; // std::vector wellperf_A_; // std::vector wellperf_phasemob_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. computeCellDynamicData(dt, state, well_state); computeFaceDynamicData(dt, state, well_state); computeWellDynamicData(dt, state, well_state); @@ -273,6 +285,8 @@ namespace Opm // std::vector cell_viscosity_; // std::vector cell_phasemob_; // std::vector cell_voldisc_; + // std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + // std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. const int nc = grid_.number_of_cells; const int np = props_.numPhases(); const double* cell_p = &state.pressure()[0]; @@ -296,6 +310,14 @@ namespace Opm // TODO: Check this! cell_voldisc_.clear(); cell_voldisc_.resize(nc, 0.0); + + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_); + rock_comp_.resize(nc); + for (int cell = 0; cell < nc; ++cell) { + rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]); + } + } } @@ -465,9 +487,16 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; - cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], - &face_gravcap_[0], cell_press, well_bhp, - &porevol_[0], h_); + if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { + cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], h_); + } else { + cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], + &face_gravcap_[0], cell_press, well_bhp, + &porevol_[0], &initial_porevol_[0], + &rock_comp_[0], h_); + } } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 8233ee17b..8a2857f21 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -33,6 +33,7 @@ namespace Opm class BlackoilState; class BlackoilPropertiesInterface; + class RockCompressibility; class LinearSolverInterface; class WellState; @@ -44,23 +45,25 @@ namespace Opm { public: /// Construct solver. - /// \param[in] grid A 2d or 3d grid. - /// \param[in] props Rock and fluid properties. - /// \param[in] linsolver Linear solver to use. - /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. - /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. - /// \param[in] maxiter Maximum acceptable number of iterations. - /// \param[in] gravity Gravity vector. If non-null, the array should - /// have D elements. - /// \param[in] wells The wells argument. Will be used in solution, - /// is ignored if NULL. - /// Note: this class observes the well object, and - /// makes the assumption that the well topology - /// and completions does not change during the - /// run. However, controls (only) are allowed - /// to change. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] props Rock and fluid properties. + /// \param[in] rock_comp_props Rock compressibility properties. May be null. + /// \param[in] linsolver Linear solver to use. + /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. + /// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller. + /// \param[in] maxiter Maximum acceptable number of iterations. + /// \param[in] gravity Gravity vector. If non-null, the array should + /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL. + /// Note: this class observes the well object, and + /// makes the assumption that the well topology + /// and completions does not change during the + /// run. However, controls (only) are allowed + /// to change. CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, const double residual_tol, const double change_tol, @@ -107,6 +110,7 @@ namespace Opm // ------ Data that will remain unmodified after construction. ------ const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; + const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; const double residual_tol_; const double change_tol_; @@ -115,7 +119,6 @@ namespace Opm const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve(). std::vector htrans_; std::vector trans_ ; - std::vector porevol_; std::vector allcells_; // ------ Internal data for the cfs_tpfa_res solver. ------ @@ -123,6 +126,7 @@ namespace Opm // ------ Data that will be modified for every solve. ------ std::vector wellperf_gpot_; + std::vector initial_porevol_; // ------ Data that will be modified for every solver iteration. ------ std::vector cell_A_; @@ -135,6 +139,8 @@ namespace Opm std::vector face_gravcap_; std::vector wellperf_A_; std::vector wellperf_phasemob_; + std::vector porevol_; // Only modified if rock_comp_props_ is non-null. + std::vector rock_comp_; // Empty unless rock_comp_props_ is non-null. // The update to be applied to the pressures (cell and bhp). std::vector pressure_increment_; diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 7a52db1e5..fb8ee8795 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1213,6 +1213,79 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , } +/* ---------------------------------------------------------------------- */ +void +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G , + double dt , + struct cfs_tpfa_res_forces *forces , + const double *zc , + struct compr_quantities_gen *cq , + const double *trans , + const double *gravcap_f, + const double *cpress , + const double *wpress , + const double *porevol , + const double *porevol0 , + const double *rock_comp, + struct cfs_tpfa_res_data *h ) +/* ---------------------------------------------------------------------- */ +{ + /* We want to add this term to the usual residual: + * + * (porevol(pressure)-porevol(initial_pressure))/dt. + * + * Its derivative (for the diagonal term of the Jacobian) is: + * + * porevol(pressure)*rock_comp(pressure)/dt + */ + + int c, w, well_is_neumann, rock_is_incomp; + size_t j; + double dpv; + const struct Wells* W; + + /* Assemble usual system (without rock compressibility). */ + cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + 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 + regular assembly, we undo it here. */ + if (well_is_neumann && h->pimpl->is_incomp) { + h->J->sa[0] /= 2; + } + + /* Add new terms to residual and Jacobian. */ + rock_is_incomp = 1; + for (c = 0; c < G->number_of_cells; c++) { + j = csrmatrix_elm_index(c, c, h->J); + + dpv = (porevol[c] - porevol0[c]); + if (dpv != 0.0 || rock_comp[c] != 0.0) { + rock_is_incomp = 0; + } + + h->J->sa[j] += porevol[c] * rock_comp[c]; + h->F[c] += dpv; + } + + /* Re-do the singularity-removing adjustment if necessary */ + if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) { + h->J->sa[0] *= 2; + } +} + + /* ---------------------------------------------------------------------- */ void cfs_tpfa_res_flux(struct UnstructuredGrid *G , diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 4fde3ba39..6a86c61bc 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -72,6 +72,22 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, const double *porevol, struct cfs_tpfa_res_data *h); +void +cfs_tpfa_res_comprock_assemble( + struct UnstructuredGrid *G, + double dt, + struct cfs_tpfa_res_forces *forces, + const double *zc, + struct compr_quantities_gen *cq, + const double *trans, + const double *gravcap_f, + const double *cpress, + const double *wpress, + const double *porevol, + const double *porevol0, + const double *rock_comp, + struct cfs_tpfa_res_data *h); + void cfs_tpfa_res_flux(struct UnstructuredGrid *G , struct cfs_tpfa_res_forces *forces , diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index 94e3d6383..0e69bc491 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G , assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok); /* We want to solve a Newton step for the residual - * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible + * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible * */ diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index f55f67edb..06708da70 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -75,15 +75,15 @@ namespace Opm void TransportModelCompressibleTwophase::solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { darcyflux_ = darcyflux; - surfacevol0_ = surfacevol0; + surfacevol0_ = &surfacevol[0]; porevolume0_ = porevolume0; porevolume_ = porevolume; source_ = source; @@ -107,6 +107,15 @@ namespace Opm &ia_downw_[0], &ja_downw_[0]); reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + surfacevol = saturation; + const int np = props_.numPhases(); + for (int cell = 0; cell < grid_.number_of_cells; ++cell) { + for (int phase = 0; phase < np; ++phase) { + surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase]; + } + } } // Residual function r(s) for a single-cell implicit Euler transport @@ -381,6 +390,7 @@ namespace Opm std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); const int nf = grid_.number_of_faces; trans_.resize(nf); + gravflux_.resize(nf); tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); tpfa_trans_compute(const_cast(&grid_), &htrans[0], &trans_[0]); } diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 19c44b22d..542363cd3 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -54,12 +54,12 @@ namespace Opm /// \param[in, out] saturation Phase saturations. void solve(const double* darcyflux, const double* pressure, - const double* surfacevol0, const double* porevolume0, const double* porevolume, const double* source, const double dt, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); /// Initialise quantities needed by gravity solver. void initGravity();