From 547efbe7d14e7d83b4aa4436637e8d4b82b31310 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:36:57 +0200 Subject: [PATCH 01/21] Added function cfs_tpfa_res_comprock_assemble(). This is intended to handle cases with both fluid and rock compressibility. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 73 ++++++++++++++++++++++ opm/core/pressure/tpfa/cfs_tpfa_residual.h | 16 +++++ 2 files changed, 89 insertions(+) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index 7a52db1e5..ab79dbce7 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 dpvdt; + const struct Wells* W; + + /* Assemble usual system (without rock compressibility). */ + cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + cpress, wpress, porevol, 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); + + dpvdt = (porevol[c] - porevol0[c]) / dt; + if (dpvdt != 0.0 || rock_comp[c] != 0.0) { + rock_is_incomp = 0; + } + + h->J->sa[j] += porevol[c] * rock_comp[c] / dt; + h->F[c] -= dpvdt; + } + + /* 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 , From d67f49bac807b29f5e5a1470d2b8b7ea0b9305a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:38:41 +0200 Subject: [PATCH 02/21] Class CompressibleTpfa now handles rock compressibility. --- opm/core/pressure/CompressibleTpfa.cpp | 39 ++++++++++++++++++++++---- opm/core/pressure/CompressibleTpfa.hpp | 38 ++++++++++++++----------- 2 files changed, 56 insertions(+), 21 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 154f06336..dbd7ea47a 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) { + 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_; From 61fdf4a8b6fda787ab3e752082c33aa04f1c3cb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 17 Aug 2012 10:56:27 +0200 Subject: [PATCH 03/21] Bugfix: call *comprock* method only if active rock compressibility. --- opm/core/pressure/CompressibleTpfa.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index dbd7ea47a..d28db5a33 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -487,7 +487,7 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; - if (rock_comp_props_ == NULL) { + 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_); From b470ab8d50e7ff8a7fc41231a86447172b3f62bf Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 20 Aug 2012 15:18:24 +0200 Subject: [PATCH 04/21] Fix initialization bug. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index f55f67edb..531aca4f5 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -381,6 +381,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]); } From 7e1e1e795138e58efea48f75e20ae7ee92dc2175 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 10:20:59 +0200 Subject: [PATCH 05/21] Fixed sign error and replaced porevol->porevol0. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index ab79dbce7..b3ca1f9d8 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1247,7 +1247,7 @@ cfs_tpfa_res_comprock_assemble( /* Assemble usual system (without rock compressibility). */ cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, - cpress, wpress, porevol, h); + cpress, wpress, porevol0, h); /* Check if we have only Neumann wells. */ well_is_neumann = 1; @@ -1276,7 +1276,7 @@ cfs_tpfa_res_comprock_assemble( } h->J->sa[j] += porevol[c] * rock_comp[c] / dt; - h->F[c] -= dpvdt; + h->F[c] += dpvdt; } /* Re-do the singularity-removing adjustment if necessary */ From 712ccb0309beafca969b66848c4ba9552b31e5e5 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:09:19 +0200 Subject: [PATCH 06/21] Corrected typo in comment. --- opm/core/pressure/tpfa/ifs_tpfa.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 * */ From 52ab67201d80a33d92f74fe9658c6792d62b9518 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:10:49 +0200 Subject: [PATCH 07/21] Changed way to get total number of cells. --- opm/core/fluid/RockFromDeck.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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); From 660f14df3b9d6cdbaf8456e273cee6a55ef7db81 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 21 Aug 2012 17:12:07 +0200 Subject: [PATCH 08/21] Corrected bug (residual should not be divided by dt). --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index b3ca1f9d8..fb8ee8795 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1242,7 +1242,7 @@ cfs_tpfa_res_comprock_assemble( int c, w, well_is_neumann, rock_is_incomp; size_t j; - double dpvdt; + double dpv; const struct Wells* W; /* Assemble usual system (without rock compressibility). */ @@ -1270,13 +1270,13 @@ cfs_tpfa_res_comprock_assemble( for (c = 0; c < G->number_of_cells; c++) { j = csrmatrix_elm_index(c, c, h->J); - dpvdt = (porevol[c] - porevol0[c]) / dt; - if (dpvdt != 0.0 || rock_comp[c] != 0.0) { + 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] / dt; - h->F[c] += dpvdt; + h->J->sa[j] += porevol[c] * rock_comp[c]; + h->F[c] += dpv; } /* Re-do the singularity-removing adjustment if necessary */ From 9ed0923279fbd33355c6ba95cee5edef148cfdb6 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 22 Aug 2012 10:11:32 +0200 Subject: [PATCH 09/21] Update surface volume in transport solver. --- .../TransportModelCompressibleTwophase.cpp | 15 ++++++++++++--- .../TransportModelCompressibleTwophase.hpp | 4 ++-- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 531aca4f5..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 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(); From 010715ad03a5c0ee67cd5be44e3ef130214175b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 Aug 2012 12:31:59 +0200 Subject: [PATCH 10/21] Document and check (primitive) for non-miscibility requirement. --- .../transport/reorder/TransportModelCompressibleTwophase.cpp | 5 +++++ .../transport/reorder/TransportModelCompressibleTwophase.hpp | 2 ++ 2 files changed, 7 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 06708da70..98aeefdf5 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -93,6 +93,11 @@ namespace Opm props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); + // Check non-miscibility requirement (only done for first cell). + if (A_[1] != 0.0 || A_[2] != 0.0) { + THROW("TransportModelCompressibleTwophase requires a property object without miscibility."); + } + std::vector seq(grid_.number_of_cells); std::vector comp(grid_.number_of_cells + 1); int ncomp; diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 542363cd3..5b82ed624 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -30,6 +30,8 @@ namespace Opm class BlackoilPropertiesInterface; + /// Implements a reordering transport solver for compressible, + /// non-miscible two-phase flow. class TransportModelCompressibleTwophase : public TransportModelInterface { public: From 96327164ccfc20db032358d172c58b805b437986 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 08:57:48 +0200 Subject: [PATCH 11/21] Added computeSurfacevol() method. The method is called by the reordering transport solver after computing new saturations in order to update the surface volumes. --- .../TransportModelCompressibleTwophase.cpp | 11 ++----- opm/core/utility/miscUtilitiesBlackoil.cpp | 30 +++++++++++++++++++ opm/core/utility/miscUtilitiesBlackoil.hpp | 19 ++++++++++++ 3 files changed, 52 insertions(+), 8 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index 98aeefdf5..b81a7e60a 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -112,15 +113,9 @@ 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]; - } - } + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } // Residual function r(s) for a single-cell implicit Euler transport diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index af088e059..1c8bceee2 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -220,5 +220,35 @@ namespace Opm } } + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases^2, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol) + { + // Note: since this is a simple matrix-vector product, it can + // be done by a BLAS call, but then we have to reorder the A + // matrix data. + std::fill(surfacevol, surfacevol + n*np, 0.0); + for (int i = 0; i < n; ++i) { + for (int row = 0; row < np; ++row) { + for (int col = 0; col < np; ++col) { + surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col]; + } + } + } + } } // namespace Opm diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 929d42ecd..4a3d6448a 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -112,6 +112,25 @@ namespace Opm const std::vector& s, std::vector& fractional_flows); + + /// Computes the surface volume densities from saturations by the formula + /// z = A s + /// for a number of data points, where z is the surface volume density, + /// s is the saturation (both as column vectors) and A is the + /// phase-to-component relation matrix. + /// @param[in] n number of data points + /// @param[in] np number of phases, must be 2 or 3 + /// @param[in] A array containing n square matrices of size num_phases, + /// in Fortran ordering, typically the output of a call + /// to the matrix() method of a BlackoilProperties* class. + /// @param[in] saturation concatenated saturation values (for all P phases) + /// @param[out] surfacevol concatenated surface-volume values (for all P phases) + void computeSurfacevol(const int n, + const int np, + const double* A, + const double* saturation, + double* surfacevol); + } // namespace Opm #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED From e96421dbd776a60c481cfda061cd8e5956691ac2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 08:59:09 +0200 Subject: [PATCH 12/21] Whitespace cleanup. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 82 +++++++++++----------- opm/core/utility/miscUtilitiesBlackoil.hpp | 26 +++---- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 1c8bceee2..548aaa23a 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -48,39 +48,39 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& press, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); std::vector visc(np); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -93,11 +93,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& press, const std::vector& z, - const std::vector& s, - std::vector& totmob) + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -126,12 +126,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -185,10 +185,10 @@ namespace Opm props.relperm(nc, &s[0], &cells[0], &pmobc[0], dpmobc); - std::transform(pmobc.begin(), pmobc.end(), - mu.begin(), - pmobc.begin(), - std::divides()); + std::transform(pmobc.begin(), pmobc.end(), + mu.begin(), + pmobc.begin(), + std::divides()); } /// Computes the fractional flow for each cell in the cells argument diff --git a/opm/core/utility/miscUtilitiesBlackoil.hpp b/opm/core/utility/miscUtilitiesBlackoil.hpp index 4a3d6448a..565acc6ef 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.hpp +++ b/opm/core/utility/miscUtilitiesBlackoil.hpp @@ -46,11 +46,11 @@ namespace Opm void computeInjectedProduced(const BlackoilPropertiesInterface& props, const std::vector& p, const std::vector& z, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced); + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced); /// @brief Computes total mobility for a set of saturation values. /// @param[in] props rock and fluid properties @@ -60,11 +60,11 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob); + const std::vector& s, + std::vector& totmob); /// @brief Computes total mobility and omega for a set of saturation values. /// @param[in] props rock and fluid properties @@ -75,12 +75,12 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props, - const std::vector& cells, + const std::vector& cells, const std::vector& p, const std::vector& z, - const std::vector& s, - std::vector& totmob, - std::vector& omega); + const std::vector& s, + std::vector& totmob, + std::vector& omega); /// @brief Computes phase mobilities for a set of saturation values. @@ -96,7 +96,7 @@ namespace Opm const std::vector& z, const std::vector& s, std::vector& pmobc); - + /// Computes the fractional flow for each cell in the cells argument /// @param[in] props rock and fluid properties From 983a29691c856444baf3b1dde998de1da1bb2cc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:00:04 +0200 Subject: [PATCH 13/21] 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). --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 42 ++++++++++------------ opm/core/pressure/tpfa/cfs_tpfa_residual.h | 13 +++++-- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index fb8ee8795..eb8dc8f34 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , /* ---------------------------------------------------------------------- */ -void +int cfs_tpfa_res_assemble(struct UnstructuredGrid *G , double dt , struct cfs_tpfa_res_forces *forces , @@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , 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); vector_zero (h->J->m, h->F); @@ -1207,14 +1207,17 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G , assemble_sources(dt, forces->src, h); } - if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] *= 2; + singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp; + if (singular) { + h->J->sa[0] *= 2.0; } + + return singular; } /* ---------------------------------------------------------------------- */ -void +int cfs_tpfa_res_comprock_assemble( struct UnstructuredGrid *G , double dt , @@ -1240,29 +1243,18 @@ cfs_tpfa_res_comprock_assemble( * porevol(pressure)*rock_comp(pressure)/dt */ - int c, w, well_is_neumann, rock_is_incomp; + int c, rock_is_incomp, singular; size_t j; - double dpv; - const struct Wells* W; + double dpv; /* 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; - } - } + singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f, + cpress, wpress, porevol0, h); /* 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; + if (singular) { + h->J->sa[0] /= 2.0; } /* Add new terms to residual and Jacobian. */ @@ -1280,9 +1272,11 @@ cfs_tpfa_res_comprock_assemble( } /* Re-do the singularity-removing adjustment if necessary */ - if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) { - h->J->sa[0] *= 2; + if (rock_is_incomp && singular) { + h->J->sa[0] *= 2.0; } + + return rock_is_incomp && singular; } diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.h b/opm/core/pressure/tpfa/cfs_tpfa_residual.h index 6a86c61bc..5eddeed68 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.h +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.h @@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G , void 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, double dt, struct cfs_tpfa_res_forces *forces, @@ -72,7 +76,12 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G, const double *porevol, 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( struct UnstructuredGrid *G, double dt, From 73949f892e91e4d07cdf63c966c593f9ca64a347 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:03:20 +0200 Subject: [PATCH 14/21] New singularPressure() method propagates singularity information. --- opm/core/pressure/CompressibleTpfa.cpp | 36 ++++++++++++++++++++------ opm/core/pressure/CompressibleTpfa.hpp | 18 +++++++++---- 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index d28db5a33..67d3a651a 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -77,7 +77,8 @@ namespace Opm wells_(wells), htrans_(grid.cell_facepos[ grid.number_of_cells ]), trans_ (grid.number_of_faces), - allcells_(grid.number_of_cells) + allcells_(grid.number_of_cells), + singular_(false) { if (wells_ && (wells_->number_of_phases != props.numPhases())) { THROW("Inconsistent number of phases specified (wells vs. props): " @@ -189,6 +190,21 @@ namespace Opm + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool CompressibleTpfa::singularPressure() const + { + return singular_; + } + + + + + /// Compute well potentials. void CompressibleTpfa::computeWellPotentials(const BlackoilState& state) { @@ -487,16 +503,20 @@ namespace Opm cq.Af = &face_A_[0]; cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; + int was_adjusted = 0; 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_); + was_adjusted = + 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_); + was_adjusted = + 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_); } + singular_ = (was_adjusted == 1); } diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 8a2857f21..9fc5c43fd 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -81,6 +81,14 @@ namespace Opm BlackoilState& state, WellState& well_state); + /// @brief After solve(), was the resulting pressure singular. + /// Returns true if the pressure is singular in the following + /// sense: if everything is incompressible and there are no + /// pressure conditions, the absolute values of the pressure + /// solution are arbitrary. (But the differences in pressure + /// are significant.) + bool singularPressure() const; + private: void computePerSolveDynamicData(const double dt, const BlackoilState& state, @@ -143,11 +151,11 @@ namespace Opm 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_; - - - - - + // True if the matrix assembled would be singular but for the + // adjustment made in the cfs_*_assemble() calls. This happens + // if everything is incompressible and there are no pressure + // conditions. + bool singular_; }; } // namespace Opm From caff665c1026d8463727c1c7b510ffd2cf4aebff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2012 14:45:23 +0200 Subject: [PATCH 15/21] Fixed solveGravity(): now properly modifies surfacevolume. Also: - solveGravity() interface changed to take surface volume as a parameter, - gravity vector is now given in initGravity() instead of solveGravity(), for consistency with the incompressible solver. --- .../TransportModelCompressibleTwophase.cpp | 22 +++++++++++++------ .../TransportModelCompressibleTwophase.hpp | 13 ++++++----- 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index b81a7e60a..c2d503275 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -39,7 +39,8 @@ namespace Opm typedef RegulaFalsi RootFinder; - TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, + TransportModelCompressibleTwophase::TransportModelCompressibleTwophase( + const UnstructuredGrid& grid, const Opm::BlackoilPropertiesInterface& props, const double tol, const int maxit) @@ -52,6 +53,7 @@ namespace Opm dt_(0.0), saturation_(grid.number_of_cells, -1.0), fractionalflow_(grid.number_of_cells, -1.0), + gravity_(0), mob_(2*grid.number_of_cells, -1.0), ia_upw_(grid.number_of_cells + 1, -1), ja_upw_(grid.number_of_faces, -1), @@ -384,7 +386,7 @@ namespace Opm - void TransportModelCompressibleTwophase::initGravity() + void TransportModelCompressibleTwophase::initGravity(const double* grav) { // Set up transmissibilities. std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); @@ -393,12 +395,15 @@ namespace Opm gravflux_.resize(nf); tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); tpfa_trans_compute(const_cast(&grid_), &htrans[0], &trans_[0]); + + // Remember gravity vector. + gravity_ = grav; } - void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav) + void TransportModelCompressibleTwophase::initGravityDynamic() { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] @@ -420,7 +425,7 @@ namespace Opm for (int ci = 0; ci < 2; ++ci) { double gdz = 0.0; for (int d = 0; d < dim; ++d) { - gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); + gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]); } gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]); } @@ -504,11 +509,11 @@ namespace Opm const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation) + std::vector& saturation, + std::vector& surfacevol) { // Assume that solve() has already been called, so that A_ is current. - initGravityDynamic(grav); + initGravityDynamic(); // Initialize mobilities. const int nc = grid_.number_of_cells; @@ -541,6 +546,9 @@ namespace Opm std::cout << "Gauss-Seidel column solver average iterations: " << double(num_iters)/double(columns.size()) << std::endl; toBothSat(saturation_, saturation); + + // Compute surface volume as a postprocessing step from saturation and A_ + computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]); } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index 5b82ed624..4ee93e0af 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -54,6 +54,7 @@ namespace Opm /// \param[in] source Transport source term. /// \param[in] dt Time step. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solve(const double* darcyflux, const double* pressure, const double* porevolume0, @@ -64,7 +65,8 @@ namespace Opm std::vector& surfacevol); /// Initialise quantities needed by gravity solver. - void initGravity(); + /// \param[in] grav Gravity vector + void initGravity(const double* grav); /// Solve for gravity segregation. /// This uses a column-wise nonlinear Gauss-Seidel approach. @@ -74,14 +76,14 @@ namespace Opm /// \param[in] columns Vector of cell-columns. /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] dt Time step. - /// \param[in] grav Gravity vector. /// \param[in, out] saturation Phase saturations. + /// \param[in, out] surfacevol Surface volume densities for each phase. void solveGravity(const std::vector >& columns, const double* pressure, const double* porevolume0, const double dt, - const double* grav, - std::vector& saturation); + std::vector& saturation, + std::vector& surfacevol); private: virtual void solveSingleCell(const int cell); @@ -90,7 +92,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - void initGravityDynamic(const double* grav); + void initGravityDynamic(); private: const UnstructuredGrid& grid_; @@ -112,6 +114,7 @@ namespace Opm std::vector saturation_; // P (= num. phases) per cell std::vector fractionalflow_; // = m[0]/(m[0] + m[1]) per cell // For gravity segregation. + const double* gravity_; std::vector trans_; std::vector density_; std::vector gravflux_; From c7bbf1146cb126e45e05b122601d0601d130320e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 08:18:59 +0200 Subject: [PATCH 16/21] Removed unneeded function numGlobalCells(). --- opm/core/fluid/RockFromDeck.cpp | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 02f675567..7ba8903fd 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -36,8 +36,6 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::tr1::array& kmap); - - int numGlobalCells(const EclipseGridParser& parser); } // anonymous namespace @@ -334,26 +332,6 @@ namespace Opm return kind; } - int numGlobalCells(const EclipseGridParser& parser) - { - int ngc = -1; - - if (parser.hasField("DIMENS")) { - const std::vector& - dims = parser.getIntegerValue("DIMENS"); - - ngc = dims[0] * dims[1] * dims[2]; - } - else if (parser.hasField("SPECGRID")) { - const SPECGRID& sgr = parser.getSPECGRID(); - - ngc = sgr.dimensions[ 0 ]; - ngc *= sgr.dimensions[ 1 ]; - ngc *= sgr.dimensions[ 2 ]; - } - - return ngc; - } } // anonymous namespace } // namespace Opm From 46fb4884106883fb5e8d14d26c0531290bf8279a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 13:31:23 +0200 Subject: [PATCH 17/21] Minor code cleanup in TransportModelTwophase. --- .../transport/reorder/TransportModelTwophase.cpp | 14 +++++++++----- .../transport/reorder/TransportModelTwophase.hpp | 9 +++++---- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b6451809b..2740d065e 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -52,8 +52,8 @@ namespace Opm source_(0), dt_(0.0), saturation_(grid.number_of_cells, -1.0), - reorder_iterations_(grid.number_of_cells, 0), fractionalflow_(grid.number_of_cells, -1.0), + reorder_iterations_(grid.number_of_cells, 0), mob_(2*grid.number_of_cells, -1.0) #ifdef EXPERIMENT_GAUSS_SEIDEL , ia_upw_(grid.number_of_cells + 1, -1), @@ -107,6 +107,13 @@ namespace Opm toBothSat(saturation_, saturation); } + + const std::vector& TransportModelTwophase::getReorderIterations() const + { + return reorder_iterations_; + } + + // Residual function r(s) for a single-cell implicit Euler transport // // r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) @@ -645,10 +652,7 @@ namespace Opm toBothSat(saturation_, saturation); } - void TransportModelTwophase::getReorderIterations() - { - return reorder_iterations_; - }; + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 358baaa31..3630ada8f 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -74,10 +74,11 @@ namespace Opm const double* porevolume, const double dt, std::vector& saturation); - //// Return reorder iterations - //// + + //// Return the number of iterations used by the reordering solver. //// \param[out] vector of iteration per cell - const std::vector& getReorderIterations(){return reorder_iterations_;}; + const std::vector& getReorderIterations() const; + private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); @@ -86,7 +87,7 @@ namespace Opm const int pos, const double* gravflux); int solveGravityColumn(const std::vector& cells); - private: + private: const UnstructuredGrid& grid_; const IncompPropertiesInterface& props_; const double* visc_; From 0da049f4c8cb1b37f91119994b8830f7d59a42a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Fri, 24 Aug 2012 13:39:42 +0200 Subject: [PATCH 18/21] Corrected mistake from moving function calls. --- opm/core/transport/reorder/TransportModelTwophase.cpp | 2 +- opm/core/transport/reorder/TransportModelTwophase.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index b6451809b..3346d28ea 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -645,7 +645,7 @@ namespace Opm toBothSat(saturation_, saturation); } - void TransportModelTwophase::getReorderIterations() + const std::vector& TransportModelTwophase::getReorderIterations() { return reorder_iterations_; }; diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 358baaa31..64449f9f1 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -77,7 +77,7 @@ namespace Opm //// Return reorder iterations //// //// \param[out] vector of iteration per cell - const std::vector& getReorderIterations(){return reorder_iterations_;}; + const std::vector& getReorderIterations(); private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); From 0a8ac1ddb52bd9e92daee9b1fc79f6c936d1c1f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 09:48:06 +0200 Subject: [PATCH 19/21] Minor revision, mostly whitespace cleanup and comments. --- opm/core/pressure/CompressibleTpfa.cpp | 6 +++--- opm/core/pressure/CompressibleTpfa.hpp | 16 ++++++++-------- .../TransportModelCompressibleTwophase.cpp | 2 +- .../transport/reorder/TransportModelTwophase.hpp | 2 +- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 67d3a651a..954e07fa4 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -504,13 +504,13 @@ namespace Opm cq.phasemobf = &face_phasemob_[0]; cq.voldiscr = &cell_voldisc_[0]; int was_adjusted = 0; - if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) { - was_adjusted = + if (! (rock_comp_props_ && rock_comp_props_->isActive())) { + was_adjusted = cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0], &face_gravcap_[0], cell_press, well_bhp, &porevol_[0], h_); } else { - was_adjusted = + was_adjusted = cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0], &face_gravcap_[0], cell_press, well_bhp, &porevol_[0], &initial_porevol_[0], diff --git a/opm/core/pressure/CompressibleTpfa.hpp b/opm/core/pressure/CompressibleTpfa.hpp index 9fc5c43fd..054ac2e4f 100644 --- a/opm/core/pressure/CompressibleTpfa.hpp +++ b/opm/core/pressure/CompressibleTpfa.hpp @@ -61,7 +61,7 @@ namespace Opm /// and completions does not change during the /// run. However, controls (only) are allowed /// to change. - CompressibleTpfa(const UnstructuredGrid& grid, + CompressibleTpfa(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, const LinearSolverInterface& linsolver, @@ -71,8 +71,8 @@ namespace Opm const double* gravity, const Wells* wells); - /// Destructor. - ~CompressibleTpfa(); + /// Destructor. + ~CompressibleTpfa(); /// Solve the pressure equation by Newton-Raphson scheme. /// May throw an exception if the number of iterations @@ -112,11 +112,11 @@ namespace Opm void solveIncrement(); double residualNorm() const; double incrementNorm() const; - void computeResults(BlackoilState& state, + void computeResults(BlackoilState& state, WellState& well_state) const; // ------ Data that will remain unmodified after construction. ------ - const UnstructuredGrid& grid_; + const UnstructuredGrid& grid_; const BlackoilPropertiesInterface& props_; const RockCompressibility* rock_comp_props_; const LinearSolverInterface& linsolver_; @@ -125,12 +125,12 @@ namespace Opm const int maxiter_; const double* gravity_; // May be NULL const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve(). - std::vector htrans_; - std::vector trans_ ; + std::vector htrans_; + std::vector trans_ ; std::vector allcells_; // ------ Internal data for the cfs_tpfa_res solver. ------ - struct cfs_tpfa_res_data* h_; + struct cfs_tpfa_res_data* h_; // ------ Data that will be modified for every solve. ------ std::vector wellperf_gpot_; diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index c2d503275..697f94a57 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -96,7 +96,7 @@ namespace Opm props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); - // Check non-miscibility requirement (only done for first cell). + // Check immiscibility requirement (only done for first cell). if (A_[1] != 0.0 || A_[2] != 0.0) { THROW("TransportModelCompressibleTwophase requires a property object without miscibility."); } diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 3630ada8f..97386c8af 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -76,7 +76,7 @@ namespace Opm std::vector& saturation); //// Return the number of iterations used by the reordering solver. - //// \param[out] vector of iteration per cell + //// \return vector of iteration per cell const std::vector& getReorderIterations() const; private: From 9bb76d74cc7987ab9696d1f135976edf6d2fc6de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 11:19:22 +0200 Subject: [PATCH 20/21] Fixed bug in matrix multiplication (matrix has Fortran element order). --- opm/core/utility/miscUtilitiesBlackoil.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 548aaa23a..74eb1c9c5 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -245,7 +245,7 @@ namespace Opm for (int i = 0; i < n; ++i) { for (int row = 0; row < np; ++row) { for (int col = 0; col < np; ++col) { - surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col]; + surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col]; } } } From 5182fef48d18ee6b47f41578c57e811dcbb83c5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 13:17:27 +0200 Subject: [PATCH 21/21] Switch loop ordering for better cache performance. --- opm/core/utility/miscUtilitiesBlackoil.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/utility/miscUtilitiesBlackoil.cpp b/opm/core/utility/miscUtilitiesBlackoil.cpp index 74eb1c9c5..6f4f75565 100644 --- a/opm/core/utility/miscUtilitiesBlackoil.cpp +++ b/opm/core/utility/miscUtilitiesBlackoil.cpp @@ -243,8 +243,8 @@ namespace Opm // matrix data. std::fill(surfacevol, surfacevol + n*np, 0.0); for (int i = 0; i < n; ++i) { - for (int row = 0; row < np; ++row) { - for (int col = 0; col < np; ++col) { + for (int col = 0; col < np; ++col) { + for (int row = 0; row < np; ++row) { surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col]; } }