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

@@ -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();

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();

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -39,7 +39,8 @@ namespace Opm
     typedef RegulaFalsi<WarnAndContinueOnError> 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<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
@@ -393,12 +395,15 @@ namespace Opm
         gravflux_.resize(nf);
         tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
         tpfa_trans_compute(const_cast<UnstructuredGrid*>(&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<double>& saturation,
-                                                          std::vector<double>& saturation)
+                                                          std::vector<double>& 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

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<double>& 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<std::vector<int> >& columns,
                           const double* pressure,
                           const double* porevolume0,
                           const double dt,
-                          const double* grav,
+                          std::vector<double>& saturation,
-                          std::vector<double>& saturation);
+                          std::vector<double>& 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<int>& cells);
-        void initGravityDynamic(const double* grav);
+        void initGravityDynamic();
 
     private:
         const UnstructuredGrid& grid_;
@@ -112,6 +114,7 @@ namespace Opm
         std::vector<double> saturation_;        // P (= num. phases) per cell
         std::vector<double> fractionalflow_;  // = m[0]/(m[0] + m[1]) per cell
         // For gravity segregation.
+        const double* gravity_;
         std::vector<double> trans_;
         std::vector<double> density_;
         std::vector<double> gravflux_;