Now running with rock compressibility (testing in progress). Multiple changes.

- TransportModelTwophase no longer takes pore volume in constructor, but in
   the solve() and solveGravity() calls.
 - Residual function uses compressibility term (not yet for gravity residual).
 - spu_2p now takes a new parameter "init_p_bar", and ReservoirState class
   accepts initial pressure as a constructor argument.
 - Moved parts of initialization around, since pore volume now depends on
   state (pressure).

examples/spu_2p.cpp

@@ -99,7 +99,7 @@
 class ReservoirState
 {
 public:
-    ReservoirState(const UnstructuredGrid* g, const double init_sat = 0.0)
+    ReservoirState(const UnstructuredGrid* g, const double init_sat, const double init_p)
         : press_ (g->number_of_cells, 0.0),
           fpress_(g->number_of_faces, 0.0),
           flux_  (g->number_of_faces, 0.0),
@@ -108,6 +108,7 @@ public:
         for (int cell = 0; cell < g->number_of_cells; ++cell) {
             sat_[2*cell] = init_sat;
             sat_[2*cell + 1] = 1.0 - init_sat;
+            press_[cell] = init_p;
         }
     }
 
@@ -425,11 +426,6 @@ main(int argc, char** argv)
         rock_comp.reset(new Opm::RockCompressibility(param));
     }
 
-    // Extra rock init.
-    std::vector<double> porevol;
-    computePorevolume(*grid->c_grid(), *props, porevol);
-    double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
-
     // Extra fluid init for transport solver.
     TwophaseFluid fluid(*props);
 
@@ -456,48 +452,39 @@ main(int argc, char** argv)
         }
     }
 
-    // Solvers init.
+    // Check that rock compressibility is not used with solvers that do not handle it.
-    // Pressure solver.
+    if (rock_comp->isActive()) {
-#ifdef EXPERIMENT_ISTL
+        if (!use_reorder) {
-    Opm::LinearSolverIstl linsolver(param);
+            THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
-#else
+        }
-    Opm::LinearSolverUmfpack linsolver;
+        if (use_segregation_split) {
-#endif // EXPERIMENT_ISTL
+            if (!use_gauss_seidel_gravity) {
-    const double *grav = use_gravity ? &gravity[0] : 0;
+                THROW("For gravity segregation splitting, only use_gauss_seidel_gravity=true supports rock compressibility.");
-    Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver);
-    // Reordering solver.
-    const double nltol = param.getDefault("nl_tolerance", 1e-9);
-    const int maxit = param.getDefault("nl_maxiter", 30);
-    Opm::TransportModelTwophase reorder_model(*grid->c_grid(), &porevol[0], *props, nltol, maxit);
-    if (use_gauss_seidel_gravity) {
-        reorder_model.initGravity(grav);
             }
-    // Non-reordering solver.
-    TransportModel  model  (fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
-    if (use_gravity) {
-        model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans());
         }
-    TransportSolver tsolver(model);
-    // Column-based gravity segregation solver.
-    typedef std::map<int, std::vector<int> > ColMap;
-    ColMap columns;
-    if (use_column_solver) {
-        Opm::extractColumn(*grid->c_grid(), columns);
     }
-    Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid(), nltol, maxit);
-
-    // Boundary conditions.
-    Opm::FlowBCManager bcs;
 
     // State-related and source-related variables init.
     int num_cells = grid->c_grid()->number_of_cells;
     std::vector<double> totmob;
     std::vector<double> omega; // Will remain empty if no gravity.
+    std::vector<double> rc; // Will remain empty if no rock compressibility.
     double init_sat = param.getDefault("init_sat", 0.0);
-    ReservoirState state(grid->c_grid(), init_sat);
+    double init_p = param.getDefault("init_p_bar", 235)*Opm::unit::barsa;
+    ReservoirState state(grid->c_grid(), init_sat, init_p);
     if (!param.has("init_sat")) {
         state.setToMinimumWaterSat(*props);
     }
+
+    // Extra rock init.
+    std::vector<double> porevol;
+    if (rock_comp->isActive()) {
+        computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
+    } else {
+        computePorevolume(*grid->c_grid(), *props, porevol);
+    }
+    double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+
     // We need a separate reorder_sat, because the reorder
     // code expects a scalar sw, not both sw and so.
     std::vector<double> reorder_sat(num_cells);
@@ -613,13 +600,44 @@ main(int argc, char** argv)
     }
     std::vector<double> reorder_src = src;
 
-    // Dirichlet boundary conditions.
+    // Boundary conditions.
+    Opm::FlowBCManager bcs;
     if (param.getDefault("use_pside", false)) {
         int pside = param.get<int>("pside");
         double pside_pressure = param.get<double>("pside_pressure");
         bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure);
     }
 
+    // Solvers init.
+    // Pressure solver.
+#ifdef EXPERIMENT_ISTL
+    Opm::LinearSolverIstl linsolver(param);
+#else
+    Opm::LinearSolverUmfpack linsolver;
+#endif // EXPERIMENT_ISTL
+    const double *grav = use_gravity ? &gravity[0] : 0;
+    Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver);
+    // Reordering solver.
+    const double nltol = param.getDefault("nl_tolerance", 1e-9);
+    const int maxit = param.getDefault("nl_maxiter", 30);
+    Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nltol, maxit);
+    if (use_gauss_seidel_gravity) {
+        reorder_model.initGravity(grav);
+    }
+    // Non-reordering solver.
+    TransportModel  model  (fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
+    if (use_gravity) {
+        model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans());
+    }
+    TransportSolver tsolver(model);
+    // Column-based gravity segregation solver.
+    typedef std::map<int, std::vector<int> > ColMap;
+    ColMap columns;
+    if (use_column_solver) {
+        Opm::extractColumn(*grid->c_grid(), columns);
+    }
+    Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid(), nltol, maxit);
+
     // Control init.
     Opm::ImplicitTransportDetails::NRReport  rpt;
     Opm::ImplicitTransportDetails::NRControl ctrl;
@@ -685,7 +703,17 @@ main(int argc, char** argv)
             computeTotalMobility(*props, allcells, state.saturation(), totmob);
         }
         pressure_timer.start();
+        if (rock_comp->isActive()) {
+            rc.resize(num_cells);
+            for (int cell = 0; cell < num_cells; ++cell) {
+                rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
+            }
+            psolver.solve(totmob, omega, src, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
+                          state.pressure(), state.faceflux());
+            computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
+        } else {
             psolver.solve(totmob, omega, src, bcs.c_bcs(), state.pressure(), state.faceflux());
+        }
         pressure_timer.stop();
         double pt = pressure_timer.secsSinceStart();
         std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
@@ -709,13 +737,14 @@ main(int argc, char** argv)
         transport_timer.start();
         if (use_reorder) {
             Opm::toWaterSat(state.saturation(), reorder_sat);
-            reorder_model.solve(&state.faceflux()[0], &reorder_src[0], simtimer.currentStepLength(), &reorder_sat[0]);
+            reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0],
+                                simtimer.currentStepLength(), &reorder_sat[0]);
             Opm::toBothSat(reorder_sat, state.saturation());
             Opm::computeInjectedProduced(*props, state.saturation(), src, simtimer.currentStepLength(), injected, produced);
             if (use_segregation_split) {
                 if (use_column_solver) {
                     if (use_gauss_seidel_gravity) {
-                        reorder_model.solveGravity(columns, simtimer.currentStepLength(), reorder_sat);
+                        reorder_model.solveGravity(columns, &porevol[0], simtimer.currentStepLength(), reorder_sat);
                         Opm::toBothSat(reorder_sat, state.saturation());
                     } else {
                         colsolver.solve(columns, simtimer.currentStepLength(), state.saturation());

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -38,12 +38,10 @@ namespace Opm
 
 
     TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid,
-						   const double* porevolume,
 						   const Opm::IncompPropertiesInterface& props,
 						   const double tol,
 						   const int maxit)
 	: grid_(grid),
-	  porevolume_(porevolume),
 	  props_(props),
 	  tol_(tol),
 	  maxit_(maxit),
@@ -75,11 +73,13 @@ namespace Opm
     }
 
     void TransportModelTwophase::solve(const double* darcyflux,
+                                       const double* porevolume,
 				       const double* source,
 				       const double dt,
 				       double* saturation)
     {
 	darcyflux_ = darcyflux;
+        porevolume_ = porevolume;
 	source_ = source;
 	dt_ = dt;
 	saturation_ = saturation;
@@ -112,8 +112,9 @@ namespace Opm
     {
 	int cell;
 	double s0;
-	double influx;  // sum_j min(v_ij, 0)*f(s_j)
+	double influx;    // sum_j min(v_ij, 0)*f(s_j) + q_w
-	double outflux; // sum_j max(v_ij, 0)
+	double outflux;   // sum_j max(v_ij, 0) - q
+        double comp_term; // q - sum_j v_ij
 	double dtpv;    // dt/pv(i)
 	const TransportModelTwophase& tm;
 	explicit Residual(const TransportModelTwophase& tmodel, int cell_index)
@@ -121,10 +122,11 @@ namespace Opm
 	{
 	    cell    = cell_index;
 	    s0      = tm.saturation_[cell];
-            double dflux       = -tm.source_[cell];
+            double src_flux       = -tm.source_[cell];
-            bool src_is_inflow = dflux < 0.0;
+            bool src_is_inflow = src_flux < 0.0;
-	    influx  =  src_is_inflow ? dflux : 0.0;
+	    influx  =  src_is_inflow ? src_flux : 0.0;
-	    outflux = !src_is_inflow ? dflux : 0.0;
+	    outflux = !src_is_inflow ? src_flux : 0.0;
+            comp_term = tm.source_[cell];   // Note: this assumes that all source flux is water.
 	    dtpv    = tm.dt_/tm.porevolume_[cell];
 
 	    for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
@@ -146,12 +148,13 @@ namespace Opm
 		    } else {
 			outflux += flux;
 		    }
+                    comp_term -= flux;
 		}
 	    }
 	}
 	double operator()(double s) const
 	{
-	    return s - s0 +  dtpv*(outflux*tm.fracFlow(s, cell) + influx);
+	    return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + dtpv*s*comp_term;
 	}
     };
 
@@ -597,6 +600,7 @@ namespace Opm
 
 
     void TransportModelTwophase::solveGravity(const std::map<int, std::vector<int> >& columns,
+                                              const double* porevolume,
                                               const double dt,
                                               std::vector<double>& saturation)
     {
@@ -617,6 +621,7 @@ namespace Opm
         }
 
         // Set up other variables.
+        porevolume_ = porevolume;
         dt_ = dt;
         saturation_ = &saturation[0];
 

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -35,12 +35,12 @@ namespace Opm
     {
     public:
 	TransportModelTwophase(const UnstructuredGrid& grid,
-			       const double* porevolume,
 			       const Opm::IncompPropertiesInterface& props,
 			       const double tol,
 			       const int maxit);
 
 	void solve(const double* darcyflux,
+                   const double* porevolume,
 		   const double* source,
 		   const double dt,
 		   double* saturation);
@@ -54,12 +54,12 @@ namespace Opm
                                     const double* gravflux);
         void solveGravityColumn(const std::vector<int>& cells);
         void solveGravity(const std::map<int, std::vector<int> >& columns,
+                          const double* porevolume,
                           const double dt,
                           std::vector<double>& saturation);
 
     private:
 	const UnstructuredGrid& grid_;
-	const double* porevolume_;  // one volume per cell
 	const IncompPropertiesInterface& props_;
 	const double* visc_;
 	std::vector<double> smin_;
@@ -68,6 +68,7 @@ namespace Opm
 	double maxit_;
 
 	const double* darcyflux_;   // one flux per grid face
+	const double* porevolume_;  // one volume per cell
 	const double* source_;      // one source per cell
 	double dt_;
 	double* saturation_;        // one per cell