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] 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_;