Class CompressibleTpfa now handles rock compressibility.

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