diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 249cbcf09..cdf443f5d 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -187,6 +188,8 @@ try Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); + std::vector threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid()); + SimulatorFullyImplicitBlackoil 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; diff --git a/examples/sim_fibo_ad_cp.cpp b/examples/sim_fibo_ad_cp.cpp index 10339668b..9b21c1457 100644 --- a/examples/sim_fibo_ad_cp.cpp +++ b/examples/sim_fibo_ad_cp.cpp @@ -41,8 +41,7 @@ #include #include - -#include +#include #include #include @@ -53,6 +52,7 @@ #include #include #include +#include // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix. #include #include @@ -68,7 +68,6 @@ #include #include -#include #include #include @@ -84,7 +83,6 @@ #include #include - namespace { void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) @@ -225,6 +223,8 @@ try Opm::DerivedGeology geology(*grid, *new_props, eclipseState, grav); + std::vector threshold_pressures = thresholdPressures(deck, eclipseState, *grid); + SimulatorFullyImplicitBlackoil 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; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index c96f898be..0fe6c2873 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/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& 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 rq_; std::vector phaseCondition_; @@ -223,6 +235,8 @@ namespace Opm { const ADB& p , const SolutionState& state ); + void applyThresholdPressures(ADB& dp); + double residualNorm() const; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 167776cd5..94a678b80 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/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(fluid.numPhases(), ADB::null()), @@ -213,6 +214,28 @@ namespace { } + + template + void + FullyImplicitBlackoilSolver:: + setThresholdPressures(const std::vector& 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 void FullyImplicitBlackoilSolver:: @@ -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 + void + FullyImplicitBlackoilSolver::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(); + + // 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(); + 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 double FullyImplicitBlackoilSolver::residualNorm() const diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 4b320918d..a33356424 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/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 eclipse_state, - EclipseWriter& output_writer); + EclipseWriter& output_writer, + const std::vector& threshold_pressures_by_face); /// Run the simulation. /// This will run succesive timesteps until timer.done() is true. It will diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 81eb54c78..97db1cfdc 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -83,7 +83,8 @@ namespace Opm bool has_disgas, bool has_vapoil, std::shared_ptr eclipse_state, - EclipseWriter& output_writer); + EclipseWriter& output_writer, + const std::vector& 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 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 eclipse_state, - EclipseWriter& output_writer) + EclipseWriter& output_writer, + const std::vector& 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 eclipse_state, - EclipseWriter& output_writer) + EclipseWriter& output_writer, + const std::vector& 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(AutoDiffGrid::numCells(grid_), 0)) + rateConverter_(props_, std::vector(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 SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer, BlackoilState& state) @@ -332,6 +341,9 @@ namespace Opm // Run a single step of the solver. solver_timer.start(); FullyImplicitBlackoilSolver 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();