Merge pull request #5 from atgeirr/master

Add rock compressibility to fully implicit solver.

examples/test_implicit_ad.cpp

@@ -101,7 +101,7 @@ main(int argc, char* argv[])
 
     Opm::LinearSolverFactory linsolver(param);
 
-    Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, *wells, linsolver);
+    Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver);
 
     Opm::BlackoilState state;
     initStateBasic(*g, props0, param, 0.0, state);

opm/autodiff/FullyImplicitBlackoilSolver.cpp

@@ -27,6 +27,7 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -187,11 +188,13 @@ namespace Opm {
     FullyImplicitBlackoilSolver(const UnstructuredGrid&         grid ,
                                 const BlackoilPropsAdInterface& fluid,
                                 const DerivedGeology&           geo  ,
+                                const RockCompressibility*      rock_comp_props,
                                 const Wells&                    wells,
                                 const LinearSolverInterface&    linsolver)
         : grid_  (grid)
         , fluid_ (fluid)
         , geo_   (geo)
+        , rock_comp_props_(rock_comp_props)
         , wells_ (wells)
         , linsolver_ (linsolver)
         , active_(activePhases(fluid.phaseUsage()))
@@ -523,12 +526,14 @@ namespace Opm {
         const std::vector<ADB>& sat   = state.saturation;
         const ADB&              rs    = state.rs;
 
+        const ADB pv_mult = poroMult(press);
+
         const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
         for (int phase = 0; phase < maxnp; ++phase) {
             if (active_[ phase ]) {
                 const int pos = pu.phase_pos[ phase ];
                 rq_[pos].b = fluidReciprocFVF(phase, press, rs, cells_);
-                rq_[pos].accum[aix] = rq_[pos].b * sat[pos];
+                rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
                 // DUMP(rq_[pos].b);
                 // DUMP(rq_[pos].accum[aix]);
             }
@@ -539,7 +544,7 @@ namespace Opm {
             const int po = pu.phase_pos[ Oil ];
             const int pg = pu.phase_pos[ Gas ];
 
-            rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix];
+            rq_[pg].accum[aix] += pv_mult * state.rs * rq_[po].accum[aix];
             // DUMP(rq_[pg].accum[aix]);
         }
     }
@@ -707,7 +712,7 @@ namespace Opm {
             const int oilpos = pu.phase_pos[Oil];
             const int gaspos = pu.phase_pos[Gas];
             const ADB rs_perf = subset(state.rs, well_cells);
-            well_rates_all += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np);
+            well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np);
             // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs);
             residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs;
         }
@@ -992,8 +997,9 @@ namespace Opm {
     {
         const int phase = canph_[ actph ];
         const ADB mu    = fluidViscosity(phase, state.pressure, state.rs, cells_);
+        const ADB tr_mult = transMult(state.pressure);
 
-        rq_[ actph ].mob = kr[ phase ] / mu;
+        rq_[ actph ].mob = tr_mult * kr[ phase ] / mu;
 
         const ADB rho   = fluidDensity(phase, state.pressure, state.rs, cells_);
         const ADB gflux = grav_ * rho;
@@ -1123,4 +1129,58 @@ namespace Opm {
     }
 
 
+
+
+
+    ADB
+    FullyImplicitBlackoilSolver::poroMult(const ADB& p) const
+    {
+        const int n = p.size();
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            V pm(n);
+            V dpm(n);
+            for (int i = 0; i < n; ++i) {
+                pm[i] = rock_comp_props_->poroMult(p.value()[i]);
+                dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
+            }
+            ADB::M dpm_diag = spdiag(dpm);
+            const int num_blocks = p.numBlocks();
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = dpm_diag * p.derivative()[block];
+            }
+            return ADB::function(pm, jacs);
+        } else {
+            return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
+        }
+    }
+
+
+
+
+
+    ADB
+    FullyImplicitBlackoilSolver::transMult(const ADB& p) const
+    {
+        const int n = p.size();
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            V tm(n);
+            V dtm(n);
+            for (int i = 0; i < n; ++i) {
+                tm[i] = rock_comp_props_->transMult(p.value()[i]);
+                dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
+            }
+            ADB::M dtm_diag = spdiag(dtm);
+            const int num_blocks = p.numBlocks();
+            std::vector<ADB::M> jacs(num_blocks);
+            for (int block = 0; block < num_blocks; ++block) {
+                jacs[block] = dtm_diag * p.derivative()[block];
+            }
+            return ADB::function(tm, jacs);
+        } else {
+            return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
+        }
+    }
+
+
 } // namespace Opm

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -30,6 +30,7 @@ struct Wells;
 namespace Opm {
 
     class DerivedGeology;
+    class RockCompressibility;
     class LinearSolverInterface;
     class BlackoilState;
     class WellState;
@@ -42,6 +43,7 @@ namespace Opm {
         FullyImplicitBlackoilSolver(const UnstructuredGrid&         grid ,
                                     const BlackoilPropsAdInterface& fluid,
                                     const DerivedGeology&           geo  ,
+                                    const RockCompressibility*      rock_comp_props,
                                     const Wells&                    wells,
                                     const LinearSolverInterface&    linsolver);
 
@@ -98,6 +100,7 @@ namespace Opm {
         const UnstructuredGrid&         grid_;
         const BlackoilPropsAdInterface& fluid_;
         const DerivedGeology&           geo_;
+        const RockCompressibility*      rock_comp_props_;
         const Wells&                    wells_;
         const LinearSolverInterface&    linsolver_;
         // For each canonical phase -> true if active
@@ -187,6 +190,12 @@ namespace Opm {
         ADB
         fluidRsMax(const ADB&              p,
                    const std::vector<int>& cells) const;
+
+        ADB
+        poroMult(const ADB& p) const;
+
+        ADB
+        transMult(const ADB& p) const;
     };
 } // namespace Opm
 

opm/autodiff/SimulatorFullyImplicitBlackoil.cpp

@@ -249,7 +249,7 @@ namespace Opm
           bcs_(bcs),
           gravity_(gravity),
           geo_(grid_, props_, gravity_),
-          solver_(grid_, props_, geo_, *wells_manager.c_wells(), linsolver)
+          solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
           /*                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
                                param.getDefault("nl_pressure_change_tolerance", 1.0),
                                param.getDefault("nl_pressure_maxiter", 10),