Merge branch 'rock_comp_tpfa'

opm/core/fluid/RockFromDeck.cpp

@@ -87,7 +87,7 @@ namespace Opm
                                           double perm_threshold)
     {
         const int dim              = 3;
-        const int num_global_cells = numGlobalCells(parser);
+        const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
         const int nc = grid.number_of_cells;
 
         ASSERT (num_global_cells > 0);

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]);
+        // 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];
+        if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
             cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                   &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;
 
@@ -46,6 +47,7 @@ namespace Opm
         /// Construct solver.
         /// \param[in] grid             A 2d or 3d grid.
         /// \param[in] props            Rock and fluid properties.
+        /// \param[in] rock_comp_props  Rock compressibility properties. May be null.
         /// \param[in] linsolver        Linear solver to use.
         /// \param[in] residual_tol     Solution accepted if inf-norm of residual is smaller.
         /// \param[in] change_tol       Solution accepted if inf-norm of change in pressure is smaller.
@@ -61,6 +63,7 @@ namespace Opm
         ///                                   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_;
 

opm/core/pressure/tpfa/cfs_tpfa_residual.c

@@ -1213,6 +1213,79 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
 }
 
 
+/* ---------------------------------------------------------------------- */
+void
+cfs_tpfa_res_comprock_assemble(
+                      struct UnstructuredGrid     *G        ,
+                      double                       dt       ,
+                      struct cfs_tpfa_res_forces  *forces   ,
+                      const double                *zc       ,
+                      struct compr_quantities_gen *cq       ,
+                      const double                *trans    ,
+                      const double                *gravcap_f,
+                      const double                *cpress   ,
+                      const double                *wpress   ,
+                      const double                *porevol  ,
+                      const double                *porevol0 ,
+                      const double                *rock_comp,
+                      struct cfs_tpfa_res_data    *h        )
+/* ---------------------------------------------------------------------- */
+{
+    /* We want to add this term to the usual residual:
+     *
+     * (porevol(pressure)-porevol(initial_pressure))/dt.
+     *
+     * Its derivative (for the diagonal term of the Jacobian) is:
+     *
+     * porevol(pressure)*rock_comp(pressure)/dt
+     */
+
+    int     c, w, well_is_neumann, rock_is_incomp;
+    size_t  j;
+    double dpv;
+    const struct Wells* W;
+
+    /* Assemble usual system (without rock compressibility). */
+    cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
+                          cpress, wpress, porevol0, h);
+
+    /* Check if we have only Neumann wells. */
+    well_is_neumann = 1;
+    W = forces->wells->W;
+    for (w = 0; well_is_neumann && w < W->number_of_wells; w++) {
+        if ((W->ctrls[w]->current >= 0) && /* OPEN? */
+            (W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) {
+            well_is_neumann = 0;
+        }
+    }
+
+    /* If we made a singularity-removing adjustment in the
+       regular assembly, we undo it here. */
+    if (well_is_neumann && h->pimpl->is_incomp) {
+        h->J->sa[0] /= 2;
+    }
+
+    /* Add new terms to residual and Jacobian. */
+    rock_is_incomp = 1;
+    for (c = 0; c < G->number_of_cells; c++) {
+        j = csrmatrix_elm_index(c, c, h->J);
+
+        dpv = (porevol[c] - porevol0[c]);
+        if (dpv != 0.0 || rock_comp[c] != 0.0) {
+            rock_is_incomp = 0;
+        }
+
+        h->J->sa[j] += porevol[c] * rock_comp[c];
+        h->F[c]     += dpv;
+    }
+
+    /* Re-do the singularity-removing adjustment if necessary */
+    if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) {
+        h->J->sa[0] *= 2;
+    }
+}
+
+
 /* ---------------------------------------------------------------------- */
 void
 cfs_tpfa_res_flux(struct UnstructuredGrid    *G        ,

opm/core/pressure/tpfa/cfs_tpfa_residual.h

@@ -72,6 +72,22 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G,
                       const double                *porevol,
                       struct cfs_tpfa_res_data    *h);
 
+void
+cfs_tpfa_res_comprock_assemble(
+                      struct UnstructuredGrid     *G,
+                      double                       dt,
+                      struct cfs_tpfa_res_forces  *forces,
+                      const double                *zc,
+                      struct compr_quantities_gen *cq,
+                      const double                *trans,
+                      const double                *gravcap_f,
+                      const double                *cpress,
+                      const double                *wpress,
+                      const double                *porevol,
+                      const double                *porevol0,
+                      const double                *rock_comp,
+                      struct cfs_tpfa_res_data    *h);
+
 void
 cfs_tpfa_res_flux(struct UnstructuredGrid    *G        ,
                   struct cfs_tpfa_res_forces *forces   ,

opm/core/pressure/tpfa/ifs_tpfa.c

@@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
     assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
 
     /* We want to solve a Newton step for the residual
-     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible
+     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible
      *
      */
 

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -75,15 +75,15 @@ namespace Opm
 
     void TransportModelCompressibleTwophase::solve(const double* darcyflux,
                                                    const double* pressure,
-                                                   const double* surfacevol0,
                                                    const double* porevolume0,
                                                    const double* porevolume,
                                                    const double* source,
                                                    const double dt,
-                                                   std::vector<double>& saturation)
+                                                   std::vector<double>& saturation,
+                                                   std::vector<double>& surfacevol)
     {
         darcyflux_ = darcyflux;
-        surfacevol0_ = surfacevol0;
+        surfacevol0_ = &surfacevol[0];
         porevolume0_ = porevolume0;
         porevolume_ = porevolume;
         source_ = source;
@@ -107,6 +107,15 @@ namespace Opm
                                &ia_downw_[0], &ja_downw_[0]);
         reorderAndTransport(grid_, darcyflux);
         toBothSat(saturation_, saturation);
+        
+        // Compute surface volume as a postprocessing step from saturation and A_
+        surfacevol = saturation;
+        const int np = props_.numPhases();
+        for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
+            for (int phase = 0; phase < np; ++phase) {
+                surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase];
+            }
+        } 
     }
 
     // Residual function r(s) for a single-cell implicit Euler transport
@@ -381,6 +390,7 @@ namespace Opm
         std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
         const int nf = grid_.number_of_faces;
         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], &trans_[0]);
     }

opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp

@@ -54,12 +54,12 @@ namespace Opm
         /// \param[in, out] saturation   Phase saturations.
         void solve(const double* darcyflux,
                    const double* pressure,
-                   const double* surfacevol0,
                    const double* porevolume0,
                    const double* porevolume,
                    const double* source,
                    const double dt,
-                   std::vector<double>& saturation);
+                   std::vector<double>& saturation,
+                   std::vector<double>& surfacevol);
 
         /// Initialise quantities needed by gravity solver.
         void initGravity();