Merge pull request #190 from atgeirr/threshold-pressure

Threshold pressure

examples/sim_fibo_ad.cpp

@@ -31,6 +31,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/thresholdPressures.hpp>
 
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/props/BlackoilPropertiesBasic.hpp>
@@ -187,6 +188,8 @@ try
 
     Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
 
+    std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid());
+
     SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
                                              *grid->c_grid(),
                                              geology,
@@ -197,7 +200,8 @@ try
                                              deck->hasKeyword("DISGAS"),
                                              deck->hasKeyword("VAPOIL"),
                                              eclipseState,
-                                             outputWriter);
+                                             outputWriter,
+                                             threshold_pressures);
 
     std::cout << "\n\n================ Starting main simulation loop ===============\n"
               << std::flush;

examples/sim_fibo_ad_cp.cpp

@@ -41,8 +41,7 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/grid/cornerpoint_grid.h>
-
+#include <opm/autodiff/GridHelpers.hpp>
-#include <opm/core/grid/GridManager.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/WellsManager.hpp>
@@ -53,6 +52,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
 
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/props/BlackoilPropertiesBasic.hpp>
@@ -68,7 +68,6 @@
 
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
-#include <opm/autodiff/GridHelpers.hpp>
 #include <opm/core/utility/share_obj.hpp>
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -84,7 +83,6 @@
 #include <vector>
 #include <numeric>
 
-
 namespace
 {
     void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
@@ -225,6 +223,8 @@ try
 
     Opm::DerivedGeology geology(*grid, *new_props, eclipseState, grav);
 
+    std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid);
+
     SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
                                                            *grid,
                                                            geology,
@@ -235,7 +235,8 @@ try
                                                            deck->hasKeyword("DISGAS"),
                                                            deck->hasKeyword("VAPOIL"),
                                                            eclipseState,
-                                                           outputWriter);
+                                                           outputWriter,
+                                                           threshold_pressures);
 
     std::cout << "\n\n================ Starting main simulation loop ===============\n"
               << std::flush;

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -74,6 +74,16 @@ namespace Opm {
                                     const bool has_disgas,
                                     const bool has_vapoil );
 
+        /// \brief Set threshold pressures that prevent or reduce flow.
+        /// This prevents flow across faces if the potential
+        /// difference is less than the threshold. If the potential
+        /// difference is greater, the threshold value is subtracted
+        /// before calculating flow. This is treated symmetrically, so
+        /// flow is prevented or reduced in both directions equally.
+        /// \param[in]  threshold_pressures_by_face   array of size equal to the number of faces
+        ///                                   of the grid passed in the constructor.
+        void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
+
         /// Take a single forward step, modifiying
         ///   state.pressure()
         ///   state.faceflux()
@@ -155,6 +165,8 @@ namespace Opm {
         double                          relax_increment_;
         double                          relax_rel_tol_;
         int                             max_iter_;
+        bool use_threshold_pressure_;
+        V threshold_pressures_by_interior_face_;
 
         std::vector<ReservoirResidualQuant> rq_;
         std::vector<PhasePresence> phaseCondition_;
@@ -223,6 +235,8 @@ namespace Opm {
                         const ADB&              p     ,
                         const SolutionState&    state );
 
+        void applyThresholdPressures(ADB& dp);
+
         double
         residualNorm() const;
 

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -190,6 +190,7 @@ namespace {
         , relax_increment_ (0.1)
         , relax_rel_tol_ (0.2)
         , max_iter_ (15)
+        , use_threshold_pressure_(false)
         , rq_    (fluid.numPhases())
         , phaseCondition_(AutoDiffGrid::numCells(grid))
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
@@ -213,6 +214,28 @@ namespace {
     }
 
 
+
+    template<class T>
+    void
+    FullyImplicitBlackoilSolver<T>::
+    setThresholdPressures(const std::vector<double>& threshold_pressures_by_face)
+    {
+        const int num_faces = AutoDiffGrid::numFaces(grid_);
+        if (int(threshold_pressures_by_face.size()) != num_faces) {
+            OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces.");
+        }
+        use_threshold_pressure_ = true;
+        // Map to interior faces.
+        const int num_ifaces = ops_.internal_faces.size();
+        threshold_pressures_by_interior_face_.resize(num_ifaces);
+        for (int ii = 0; ii < num_ifaces; ++ii) {
+            threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]];
+        }
+    }
+
+
+
+
     template<class T>
     void
     FullyImplicitBlackoilSolver<T>::
@@ -1577,7 +1600,11 @@ namespace {
         // compute gravity potensial using the face average as in eclipse and MRST
         const ADB rhoavg = ops_.caver * rho;
 
-        const ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
+        ADB dp = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
+
+        if (use_threshold_pressure_) {
+            applyThresholdPressures(dp);
+        }
 
         head = transi*dp;
         //head      = transi*(ops_.ngrad * phasePressure) + gflux;
@@ -1595,6 +1622,39 @@ namespace {
 
 
 
+    template<class T>
+    void
+    FullyImplicitBlackoilSolver<T>::applyThresholdPressures(ADB& dp)
+    {
+        // We support reversible threshold pressures only.
+        // Method: if the potential difference is lower (in absolute
+        // value) than the threshold for any face, then the potential
+        // (and derivatives) is set to zero. If it is above the
+        // threshold, the threshold pressure is subtracted from the
+        // absolute potential (the potential is moved towards zero).
+
+        // Identify the set of faces where the potential is under the
+        // threshold, that shall have zero flow. Storing the bool
+        // Array as a V (a double Array) with 1 and 0 elements, a
+        // 1 where flow is allowed, a 0 where it is not.
+        const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast<double>();
+
+        // Create a sparse vector that nullifies the low potential elements.
+        const M keep_high_potential = spdiag(high_potential);
+
+        // The threshold modification must have the sign that reduces the
+        // absolute value of the potential.
+        const V sign_for_mod = (dp.value() < 0).template cast<double>();
+        const V threshold_modification = sign_for_mod * threshold_pressures_by_interior_face_;
+
+        // Modify potential and nullify where appropriate.
+        dp = keep_high_potential * (dp + threshold_modification);
+    }
+
+
+
+
+
     template<class T>
     double
     FullyImplicitBlackoilSolver<T>::residualNorm() const

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -74,6 +74,7 @@ namespace Opm
         /// \param[in] vapoil        true for vaporized oil option
         /// \param[in] eclipse_state
         /// \param[in] output_writer
+        /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
         SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
                                        const Grid& grid,
                                        const DerivedGeology& geo,
@@ -84,7 +85,8 @@ namespace Opm
                                        const bool disgas,
                                        const bool vapoil,
                                        std::shared_ptr<EclipseState> eclipse_state,
-                                       EclipseWriter& output_writer);
+                                       EclipseWriter& output_writer,
+                                       const std::vector<double>& threshold_pressures_by_face);
 
         /// Run the simulation.
         /// This will run succesive timesteps until timer.done() is true. It will

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -83,7 +83,8 @@ namespace Opm
              bool has_disgas,
              bool has_vapoil,
              std::shared_ptr<EclipseState> eclipse_state,
-             EclipseWriter& output_writer);
+             EclipseWriter& output_writer,
+             const std::vector<double>& threshold_pressures_by_face);
 
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state);
@@ -118,6 +119,8 @@ namespace Opm
         // output_writer
         EclipseWriter& output_writer_;
         RateConverterType rateConverter_;
+        // Threshold pressures.
+        std::vector<double> threshold_pressures_by_face_;
 
         void
         computeRESV(const std::size_t               step,
@@ -140,11 +143,12 @@ namespace Opm
                                                                    const bool has_disgas,
                                                                    const bool has_vapoil,
                                                                    std::shared_ptr<EclipseState> eclipse_state,
-                                                                   EclipseWriter& output_writer)
+                                                                   EclipseWriter& output_writer,
+                                                                   const std::vector<double>& threshold_pressures_by_face)
 
     {
         pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil,
-                              eclipse_state, output_writer));
+                              eclipse_state, output_writer, threshold_pressures_by_face));
     }
 
 
@@ -229,7 +233,8 @@ namespace Opm
                                                const bool has_disgas,
                                                const bool has_vapoil,
                                                std::shared_ptr<EclipseState> eclipse_state,
-                                               EclipseWriter& output_writer)
+                                               EclipseWriter& output_writer,
+                                               const std::vector<double>& threshold_pressures_by_face)
         : param_(param),
           grid_(grid),
           props_(props),
@@ -241,7 +246,8 @@ namespace Opm
           has_vapoil_(has_vapoil),
           eclipse_state_(eclipse_state),
           output_writer_(output_writer),
-          rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0))
+          rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
+          threshold_pressures_by_face_(threshold_pressures_by_face)
     {
         // For output.
         output_ = param.getDefault("output", true);
@@ -267,6 +273,9 @@ namespace Opm
         }
     }
 
+
+
+
     template<class T>
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
                                                                  BlackoilState& state)
@@ -332,6 +341,9 @@ namespace Opm
             // Run a single step of the solver.
             solver_timer.start();
             FullyImplicitBlackoilSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
+            if (!threshold_pressures_by_face_.empty()) {
+                solver.setThresholdPressures(threshold_pressures_by_face_);
+            }
             solver.step(timer.currentStepLength(), state, well_state);
             solver_timer.stop();