Gravity segregation column solver for compressible case implemented.

examples/sim_wateroil.cpp

@@ -68,8 +68,6 @@
 #include <vector>
 #include <numeric>
 
-#define TRANSPORT_SOLVER_FIXED 1
-
 
 template <class State>
 static void outputState(const UnstructuredGrid& grid,
@@ -299,14 +297,13 @@ main(int argc, char** argv)
                                   nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
                                   grav, wells->c_wells());
     // Reordering solver.
-#if TRANSPORT_SOLVER_FIXED
     const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
     const int nl_maxiter = param.getDefault("nl_maxiter", 30);
     Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
     if (use_gauss_seidel_gravity) {
-        reorder_model.initGravity(grav);
+        reorder_model.initGravity();
     }
-#endif // TRANSPORT_SOLVER_FIXED
+
     // Column-based gravity segregation solver.
     std::vector<std::vector<int> > columns;
     if (use_column_solver) {
@@ -414,7 +411,6 @@ main(int argc, char** argv)
 
         // Solve transport.
         transport_timer.start();
-#if TRANSPORT_SOLVER_FIXED
         double stepsize = simtimer.currentStepLength();
         if (num_transport_substeps != 1) {
             stepsize /= double(num_transport_substeps);
@@ -427,11 +423,9 @@ main(int argc, char** argv)
                                 &porevol[0],  &porevol[0], &reorder_src[0], stepsize, state.saturation());
             // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
             if (use_segregation_split) {
-                THROW("Segregation not implemented yet.");
+                reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation());
-                // reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
             }
         }
-#endif // TRANSPORT_SOLVER_FIXED
         transport_timer.stop();
         double tt = transport_timer.secsSinceStart();
         std::cout << "Transport solver took: " << tt << " seconds." << std::endl;

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -254,6 +254,8 @@ namespace Opm
                 ++update_count;
                 const int cell = cells[i];
                 const double old_s = saturation_[cell];
+                // solveSingleCell() requires saturation_[cell]
+                // to be s0.
                 saturation_[cell] = s0[i];
                 solveSingleCell(cell);
                 const double s_change = std::fabs(saturation_[cell] - old_s);
@@ -305,8 +307,9 @@ namespace Opm
 
     // Residual function r(s) for a single-cell implicit Euler gravity segregation
     //
-    //     r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ).
+    // [[ incompressible was: r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij )  ]]
     //
+    //     r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) )
     struct TransportModelCompressibleTwophase::GravityResidual
     {
         int cell;
@@ -372,32 +375,52 @@ namespace Opm
 
 
 
-    void TransportModelCompressibleTwophase::initGravity(const double* grav)
+    void TransportModelCompressibleTwophase::initGravity()
     {
-        // Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
+        // Set up transmissibilities.
         std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
         const int nf = grid_.number_of_faces;
-        const int dim = grid_.dimensions;
+        trans_.resize(nf);
-        gravflux_.resize(nf);
         tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
-        tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &gravflux_[0]);
+        tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
+    }
 
-        const double delta_rho = 0.0;// props_.density()[0] - props_.density()[1];
+
-        THROW("TransportModelCompressibleTwophase gravity solver not done yet."); // See line above...
+
+
+    void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav)
+    {
+        // 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) ]
+        // But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density().
+        // We assume that we already have stored T_ij in trans_.
+        // We also assume that the A_ matrices are updated from an earlier call to solve().
+        const int nc = grid_.number_of_cells;
+        const int nf = grid_.number_of_faces;
+        const int np = props_.numPhases();
+        ASSERT(np == 2);
+        const int dim = grid_.dimensions;
+        density_.resize(nc*np);
+        props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
+        std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
         for (int f = 0; f < nf; ++f) {
             const int* c = &grid_.face_cells[2*f];
-            double gdz = 0.0;
+            const double signs[2] = { 1.0, -1.0 };
             if (c[0] != -1 && c[1] != -1) {
-                for (int d = 0; d < dim; ++d) {
+                for (int ci = 0; ci < 2; ++ci) {
-                    gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]);
+                    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]);
+                    }
+                    gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]);
                 }
             }
-            gravflux_[f] *= delta_rho*gdz;
         }
     }
 
 
 
+
     void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector<int>& cells,
                                                         const int pos,
                                                         const double* gravflux)
@@ -470,10 +493,13 @@ namespace Opm
     void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
                                                           const double* pressure,
                                                           const double* porevolume0,
-                                                          const double* porevolume,
                                                           const double dt,
+                                                          const double* grav,
                                                           std::vector<double>& saturation)
     {
+        // Assume that solve() has already been called, so that A_ is current.
+        initGravityDynamic(grav);
+
         // Initialize mobilities.
         const int nc = grid_.number_of_cells;
         std::vector<int> cells(nc);
@@ -493,7 +519,6 @@ namespace Opm
 
         // Set up other variables.
         porevolume0_ = porevolume0;
-        porevolume_ = porevolume;
         dt_ = dt;
         toWaterSat(saturation, saturation_);
 

opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp

@@ -48,7 +48,7 @@ namespace Opm
 	/// \param[in] pressure          Array of cell pressures
 	/// \param[in] surfacevol0       Array of surface volumes at start of timestep
 	/// \param[in] porevolume0       Array of pore volumes at start of timestep.
-	/// \param[in] porevolume        Array of pore volumes.
+	/// \param[in] porevolume        Array of pore volumes at end of timestep.
 	/// \param[in] source            Transport source term.
 	/// \param[in] dt                Time step.
 	/// \param[in, out] saturation   Phase saturations.
@@ -62,20 +62,23 @@ namespace Opm
 		   std::vector<double>& saturation);
 
         /// Initialise quantities needed by gravity solver.
-        /// \param[in] grav    gravity vector
+        void initGravity();
-        void initGravity(const double* grav);
 
         /// Solve for gravity segregation.
         /// This uses a column-wise nonlinear Gauss-Seidel approach.
         /// It assumes that the input columns contain cells in a single
         /// vertical stack, that do not interact with other columns (for
         /// gravity segregation.
-        /// \TODO: Implement this.
+	/// \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.
         void solveGravity(const std::vector<std::vector<int> >& columns,
                           const double* pressure,
                           const double* porevolume0,
-                          const double* porevolume,
                           const double dt,
+                          const double* grav,
                           std::vector<double>& saturation);
 
     private:
@@ -85,6 +88,7 @@ namespace Opm
                                     const int pos,
                                     const double* gravflux);
         int solveGravityColumn(const std::vector<int>& cells);
+        void initGravityDynamic(const double* grav);
 
     private:
 	const UnstructuredGrid& grid_;
@@ -106,6 +110,8 @@ 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.
+        std::vector<double> trans_;
+        std::vector<double> density_;
         std::vector<double> gravflux_;
         std::vector<double> mob_;
         std::vector<double> s0_;