Class CompressibleTpfa now handles rock compressibility.

examples/sim_wateroil.cpp

@@ -141,7 +141,6 @@ main(int argc, char** argv)
     std::cout << "---------------    Reading parameters     ---------------" << std::endl;
 
     // Reading various control parameters.
-    const bool use_reorder = param.getDefault("use_reorder", true);
     const bool output = param.getDefault("output", true);
     std::string output_dir;
     int output_interval = 1;
@@ -231,7 +230,7 @@ main(int argc, char** argv)
     bool use_segregation_split = false;
     bool use_column_solver = false;
     bool use_gauss_seidel_gravity = false;
-    if (use_gravity && use_reorder) {
+    if (use_gravity) {
         use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
         if (use_segregation_split) {
             use_column_solver = param.getDefault("use_column_solver", use_column_solver);
@@ -241,18 +240,6 @@ main(int argc, char** argv)
         }
     }
 
-    // Check that rock compressibility is not used with solvers that do not handle it.
-    // int nl_pressure_maxiter = 0;
-    // double nl_pressure_tolerance = 0.0;
-    if (rock_comp->isActive()) {
-        THROW("No rock compressibility in comp. pressure solver yet.");
-        if (!use_reorder) {
-            THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
-        }
-        // nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
-        // nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
-    }
-
     // Source-related variables init.
     int num_cells = grid->c_grid()->number_of_cells;
     std::vector<double> totmob;
@@ -262,11 +249,11 @@ main(int argc, char** argv)
     // Extra rock init.
     std::vector<double> porevol;
     if (rock_comp->isActive()) {
-        THROW("CompressibleTpfa solver does not handle this.");
         computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
     } else {
         computePorevolume(*grid->c_grid(), props->porosity(), porevol);
     }
+    std::vector<double> initial_porevol = porevol;
     double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
 
 
@@ -293,7 +280,7 @@ main(int argc, char** argv)
     const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0);
     const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20);
     const double *grav = use_gravity ? &gravity[0] : 0;
-    Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver,
+    Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
                                   nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
                                   grav, wells->c_wells());
     // Reordering solver.
@@ -409,6 +396,11 @@ main(int argc, char** argv)
         Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
                                     wells->c_wells(), well_state.perfRates(), reorder_src);
 
+        // Compute new porevolumes after pressure solve, if necessary.
+        if (rock_comp->isActive()) {
+            initial_porevol = porevol;
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+        }
         // Solve transport.
         transport_timer.start();
         double stepsize = simtimer.currentStepLength();
@@ -420,10 +412,11 @@ main(int argc, char** argv)
             // Note that for now we do not handle rock compressibility,
             // although the transport solver should be able to.
             reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0],
-                                &porevol[0],  &porevol[0], &reorder_src[0], stepsize, state.saturation());
+                                &porevol[0],  &initial_porevol[0], &reorder_src[0], stepsize, state.saturation());
             // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
             if (use_segregation_split) {
-                reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation());
+                reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
+                                           stepsize, grav, state.saturation());
             }
         }
         transport_timer.stop();

opm/core/pressure/CompressibleTpfa.cpp

@@ -30,6 +30,7 @@
 #include <opm/core/newwells.h>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/fluid/RockCompressibility.hpp>
 
 #include <algorithm>
 #include <cmath>
@@ -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<UnstructuredGrid*>(&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<double> face_gravcap_;
         // std::vector<double> wellperf_A_;
         // std::vector<double> wellperf_phasemob_;
+        // std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        // std::vector<double> 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<double> cell_viscosity_;
         // std::vector<double> cell_phasemob_;
         // std::vector<double> cell_voldisc_;
+        // std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        // std::vector<double> 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],
+        if (rock_comp_props_ == NULL) {
-                              &face_gravcap_[0], cell_press, well_bhp,
+            cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
-                              &porevol_[0], h_);
+                                  &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_);
+        }
     }
 
 

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] grid             A 2d or 3d grid.
-        /// \param[in] props         Rock and fluid properties.
+        /// \param[in] props            Rock and fluid properties.
-        /// \param[in] linsolver     Linear solver to use.
+        /// \param[in] rock_comp_props  Rock compressibility properties. May be null.
-        /// \param[in] residual_tol  Solution accepted if inf-norm of residual is smaller.
+        /// \param[in] linsolver        Linear solver to use.
-        /// \param[in] change_tol    Solution accepted if inf-norm of change in pressure is smaller.
+        /// \param[in] residual_tol     Solution accepted if inf-norm of residual is smaller.
-        /// \param[in] maxiter       Maximum acceptable number of iterations.
+        /// \param[in] change_tol       Solution accepted if inf-norm of change in pressure is smaller.
-        /// \param[in] gravity       Gravity vector. If non-null, the array should
+        /// \param[in] maxiter          Maximum acceptable number of iterations.
-        ///                          have D elements.
+        /// \param[in] gravity          Gravity vector. If non-null, the array should
-        /// \param[in] wells         The wells argument. Will be used in solution,
+        ///                             have D elements.
-        ///                          is ignored if NULL.
+        /// \param[in] wells            The wells argument. Will be used in solution,
-        ///                          Note: this class observes the well object, and
+        ///                             is ignored if NULL.
-        ///                                makes the assumption that the well topology
+        ///                             Note: this class observes the well object, and
-        ///                                and completions does not change during the
+        ///                                   makes the assumption that the well topology
-        ///                                run. However, controls (only) are allowed
+        ///                                   and completions does not change during the
-        ///                                to change.
+        ///                                   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<double> htrans_;
 	std::vector<double> trans_ ;
-        std::vector<double> porevol_;
         std::vector<int> 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<double> wellperf_gpot_;
+        std::vector<double> initial_porevol_;
 
         // ------ Data that will be modified for every solver iteration. ------
         std::vector<double> cell_A_;
@@ -135,6 +139,8 @@ namespace Opm
         std::vector<double> face_gravcap_;
         std::vector<double> wellperf_A_;
         std::vector<double> wellperf_phasemob_;
+        std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
         // The update to be applied to the pressures (cell and bhp).
         std::vector<double> pressure_increment_;