From 26149c30a2a88a95b56fe866ee62bdcd91b0fe05 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. --- examples/sim_wateroil.cpp | 4 ++-- .../SimulatorCompressibleTwophase.cpp | 5 ++--- .../TransportModelCompressibleTwophase.cpp | 22 +++++++++++++------ .../TransportModelCompressibleTwophase.hpp | 13 ++++++----- 4 files changed, 27 insertions(+), 17 deletions(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index c1f26233..f7e0172f 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -280,7 +280,7 @@ main(int argc, char** argv) const int nl_maxiter = param.getDefault("nl_maxiter", 30); Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); if (use_segregation_split) { - reorder_model.initGravity(); + reorder_model.initGravity(grav); } // Column-based gravity segregation solver. @@ -409,7 +409,7 @@ main(int argc, char** argv) // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0], - stepsize, grav, state.saturation()); + stepsize, state.saturation(), state.surfacevol()); } } transport_timer.stop(); diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 4c9d26c5..c5a77eb0 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -245,7 +245,6 @@ namespace Opm wells_(wells_manager.c_wells()), src_(src), bcs_(bcs), - gravity_(gravity), psolver_(grid, props, rock_comp, linsolver, param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), @@ -279,7 +278,7 @@ namespace Opm num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); use_segregation_split_ = param.getDefault("use_segregation_split", false); if (gravity != 0 && use_segregation_split_){ - tsolver_.initGravity(); + tsolver_.initGravity(gravity); extractColumn(grid_, columns_); } @@ -453,7 +452,7 @@ namespace Opm transport_src, stepsize, injected, produced); if (gravity_ != 0 && use_segregation_split_) { tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0], - stepsize, gravity_, state.saturation()); + stepsize, state.saturation(), state.surfacevol()); } } transport_timer.stop(); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index b81a7e60..c2d50327 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 5b82ed62..4ee93e0a 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_;