diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 990d77352..5be7b04b9 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -94,6 +94,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AutoDiffHelpers.hpp
opm/autodiff/AutoDiff.hpp
opm/autodiff/BackupRestore.hpp
+ opm/autodiff/BlackoilModel.hpp
+ opm/autodiff/BlackoilModel_impl.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/BlackoilPropsAdInterface.hpp
opm/autodiff/CPRPreconditioner.hpp
@@ -108,6 +110,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/NewtonIterationBlackoilCPR.hpp
opm/autodiff/NewtonIterationBlackoilInterface.hpp
opm/autodiff/NewtonIterationBlackoilSimple.hpp
+ opm/autodiff/NewtonSolver.hpp
+ opm/autodiff/NewtonSolver_impl.hpp
opm/autodiff/LinearisedBlackoilResidual.hpp
opm/autodiff/RateConverter.hpp
opm/autodiff/RedistributeDataHandles.hpp
diff --git a/opm/autodiff/BlackoilModel.hpp b/opm/autodiff/BlackoilModel.hpp
new file mode 100644
index 000000000..ff8317852
--- /dev/null
+++ b/opm/autodiff/BlackoilModel.hpp
@@ -0,0 +1,441 @@
+/*
+ Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+ Copyright 2014, 2015 Statoil ASA.
+ Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+ Copyright 2015 NTNU
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
+#define OPM_BLACKOILMODEL_HEADER_INCLUDED
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+ namespace parameter { class ParameterGroup; }
+ class DerivedGeology;
+ class RockCompressibility;
+ class NewtonIterationBlackoilInterface;
+ class BlackoilState;
+ class WellStateFullyImplicitBlackoil;
+
+
+ /// A model implementation for three-phase black oil.
+ ///
+ /// The simulator is capable of handling three-phase problems
+ /// where gas can be dissolved in oil and vice versa. It
+ /// uses an industry-standard TPFA discretization with per-phase
+ /// upwind weighting of mobilities.
+ ///
+ /// It uses automatic differentiation via the class AutoDiffBlock
+ /// to simplify assembly of the jacobian matrix.
+ template
+ class BlackoilModel
+ {
+ public:
+ // --------- Types and enums ---------
+ typedef AutoDiffBlock ADB;
+ typedef ADB::V V;
+ typedef ADB::M M;
+ typedef BlackoilState ReservoirState;
+ typedef WellStateFullyImplicitBlackoil WellState;
+
+ /// Model-specific solver parameters.
+ struct ModelParameters
+ {
+ double dp_max_rel_;
+ double ds_max_;
+ double dr_max_rel_;
+ double max_residual_allowed_;
+ double tolerance_mb_;
+ double tolerance_cnv_;
+ double tolerance_wells_;
+
+ explicit ModelParameters( const parameter::ParameterGroup& param );
+ ModelParameters();
+
+ void reset();
+ };
+
+ // --------- Public methods ---------
+
+ /// Construct the model. It will retain references to the
+ /// arguments of this functions, and they are expected to
+ /// remain in scope for the lifetime of the solver.
+ /// \param[in] param parameters
+ /// \param[in] grid grid data structure
+ /// \param[in] fluid fluid properties
+ /// \param[in] geo rock properties
+ /// \param[in] rock_comp_props if non-null, rock compressibility properties
+ /// \param[in] wells well structure
+ /// \param[in] linsolver linear solver
+ /// \param[in] has_disgas turn on dissolved gas
+ /// \param[in] has_vapoil turn on vaporized oil feature
+ /// \param[in] terminal_output request output to cout/cerr
+ BlackoilModel(const ModelParameters& param,
+ const Grid& grid ,
+ const BlackoilPropsAdInterface& fluid,
+ const DerivedGeology& geo ,
+ const RockCompressibility* rock_comp_props,
+ const Wells* wells,
+ const NewtonIterationBlackoilInterface& linsolver,
+ const bool has_disgas,
+ const bool has_vapoil,
+ const bool terminal_output);
+
+ /// \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);
+
+ /// Called once before each time step.
+ /// \param[in] dt time step size
+ /// \param[in, out] reservoir_state reservoir state variables
+ /// \param[in, out] well_state well state variables
+ void prepareStep(const double dt,
+ ReservoirState& reservoir_state,
+ WellState& well_state);
+
+ /// Called once after each time step.
+ /// In this class, this function does nothing.
+ /// \param[in] dt time step size
+ /// \param[in, out] reservoir_state reservoir state variables
+ /// \param[in, out] well_state well state variables
+ void afterStep(const double dt,
+ ReservoirState& reservoir_state,
+ WellState& well_state);
+
+ /// Assemble the residual and Jacobian of the nonlinear system.
+ /// \param[in] reservoir_state reservoir state variables
+ /// \param[in, out] well_state well state variables
+ /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
+ void assemble(const BlackoilState& reservoir_state,
+ WellStateFullyImplicitBlackoil& well_state,
+ const bool initial_assembly);
+
+ /// \brief Compute the residual norms of the mass balance for each phase,
+ /// the well flux, and the well equation.
+ /// \return a vector that contains for each phase the norm of the mass balance
+ /// and afterwards the norm of the residual of the well flux and the well equation.
+ std::vector computeResidualNorms() const;
+
+ /// The size (number of unknowns) of the nonlinear system of equations.
+ int sizeNonLinear() const;
+
+ /// Number of linear iterations used in last call to solveJacobianSystem().
+ int linearIterationsLastSolve() const;
+
+ /// Solve the Jacobian system Jx = r where J is the Jacobian and
+ /// r is the residual.
+ V solveJacobianSystem() const;
+
+ /// Apply an update to the primary variables, chopped if appropriate.
+ /// \param[in] dx updates to apply to primary variables
+ /// \param[in, out] reservoir_state reservoir state variables
+ /// \param[in, out] well_state well state variables
+ void updateState(const V& dx,
+ BlackoilState& reservoir_state,
+ WellStateFullyImplicitBlackoil& well_state);
+
+ /// Return true if output to cout is wanted.
+ bool terminalOutputEnabled() const;
+
+ /// Compute convergence based on total mass balance (tol_mb) and maximum
+ /// residual mass balance (tol_cnv).
+ /// \param[in] dt timestep length
+ /// \param[in] iteration current iteration number
+ bool getConvergence(const double dt, const int iteration);
+
+ /// The number of active phases in the model.
+ int numPhases() const;
+
+ private:
+
+ // --------- Types and enums ---------
+
+ typedef Eigen::Array DataBlock;
+
+ struct ReservoirResidualQuant {
+ ReservoirResidualQuant();
+ std::vector accum; // Accumulations
+ ADB mflux; // Mass flux (surface conditions)
+ ADB b; // Reciprocal FVF
+ ADB head; // Pressure drop across int. interfaces
+ ADB mob; // Phase mobility (per cell)
+ };
+
+ struct SolutionState {
+ SolutionState(const int np);
+ ADB pressure;
+ ADB temperature;
+ std::vector saturation;
+ ADB rs;
+ ADB rv;
+ ADB qs;
+ ADB bhp;
+ // Below are quantities stored in the state for optimization purposes.
+ std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
+ };
+
+ struct WellOps {
+ WellOps(const Wells* wells);
+ M w2p; // well -> perf (scatter)
+ M p2w; // perf -> well (gather)
+ };
+
+ enum { Water = BlackoilPropsAdInterface::Water,
+ Oil = BlackoilPropsAdInterface::Oil ,
+ Gas = BlackoilPropsAdInterface::Gas ,
+ MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
+ };
+
+ enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
+
+ // --------- Data members ---------
+
+ const Grid& grid_;
+ const BlackoilPropsAdInterface& fluid_;
+ const DerivedGeology& geo_;
+ const RockCompressibility* rock_comp_props_;
+ const Wells* wells_;
+ const NewtonIterationBlackoilInterface& linsolver_;
+ // For each canonical phase -> true if active
+ const std::vector active_;
+ // Size = # active phases. Maps active -> canonical phase indices.
+ const std::vector canph_;
+ const std::vector cells_; // All grid cells
+ HelperOps ops_;
+ const WellOps wops_;
+ const bool has_disgas_;
+ const bool has_vapoil_;
+
+ ModelParameters param_;
+ bool use_threshold_pressure_;
+ V threshold_pressures_by_interior_face_;
+
+ std::vector rq_;
+ std::vector phaseCondition_;
+ V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
+
+ LinearisedBlackoilResidual residual_;
+
+ /// \brief Whether we print something to std::cout
+ bool terminal_output_;
+
+ std::vector primalVariable_;
+ V pvdt_;
+
+ // --------- Private methods ---------
+
+ // return true if wells are available
+ bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
+ // return wells object
+ const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
+
+ SolutionState
+ constantState(const BlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw) const;
+
+ void
+ makeConstantState(SolutionState& state) const;
+
+ SolutionState
+ variableState(const BlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw) const;
+
+ void
+ computeAccum(const SolutionState& state,
+ const int aix );
+
+ void computeWellConnectionPressures(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw);
+
+ void
+ addWellControlEq(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw,
+ const V& aliveWells);
+
+ void
+ addWellEq(const SolutionState& state,
+ WellStateFullyImplicitBlackoil& xw,
+ V& aliveWells);
+
+ void updateWellControls(WellStateFullyImplicitBlackoil& xw) const;
+
+ std::vector
+ computePressures(const SolutionState& state) const;
+
+ std::vector
+ computePressures(const ADB& po,
+ const ADB& sw,
+ const ADB& so,
+ const ADB& sg) const;
+
+ V
+ computeGasPressure(const V& po,
+ const V& sw,
+ const V& so,
+ const V& sg) const;
+
+ std::vector
+ computeRelPerm(const SolutionState& state) const;
+
+ void
+ computeMassFlux(const int actph ,
+ const V& transi,
+ const ADB& kr ,
+ const ADB& p ,
+ const SolutionState& state );
+
+ void applyThresholdPressures(ADB& dp);
+
+ ADB
+ fluidViscosity(const int phase,
+ const ADB& p ,
+ const ADB& temp ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ ADB
+ fluidReciprocFVF(const int phase,
+ const ADB& p ,
+ const ADB& temp ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ ADB
+ fluidDensity(const int phase,
+ const ADB& p ,
+ const ADB& temp ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ V
+ fluidRsSat(const V& p,
+ const V& so,
+ const std::vector& cells) const;
+
+ ADB
+ fluidRsSat(const ADB& p,
+ const ADB& so,
+ const std::vector& cells) const;
+
+ V
+ fluidRvSat(const V& p,
+ const V& so,
+ const std::vector& cells) const;
+
+ ADB
+ fluidRvSat(const ADB& p,
+ const ADB& so,
+ const std::vector& cells) const;
+
+ ADB
+ poroMult(const ADB& p) const;
+
+ ADB
+ transMult(const ADB& p) const;
+
+ void
+ classifyCondition(const SolutionState& state,
+ std::vector& cond ) const;
+
+ const std::vector
+ phaseCondition() const {return phaseCondition_;}
+
+ void
+ classifyCondition(const BlackoilState& state);
+
+
+ /// update the primal variable for Sg, Rv or Rs. The Gas phase must
+ /// be active to call this method.
+ void
+ updatePrimalVariableFromState(const BlackoilState& state);
+
+ /// Update the phaseCondition_ member based on the primalVariable_ member.
+ void
+ updatePhaseCondFromPrimalVariable();
+
+ /// \brief Compute the reduction within the convergence check.
+ /// \param[in] B A matrix with MaxNumPhases columns and the same number rows
+ /// as the number of cells of the grid. B.col(i) contains the values
+ /// for phase i.
+ /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+ /// as the number of cells of the grid. tempV.col(i) contains the
+ /// values
+ /// for phase i.
+ /// \param[in] R A matrix with MaxNumPhases columns and the same number rows
+ /// as the number of cells of the grid. B.col(i) contains the values
+ /// for phase i.
+ /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+ /// of R for the phase i.
+ /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+ /// maximum of tempV for the phase i.
+ /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+ /// of B for the phase i.
+ /// \param[out] maxNormWell The maximum of the well equations for each phase.
+ /// \param[in] nc The number of cells of the local grid.
+ /// \param[in] nw The number of wells on the local grid.
+ /// \return The total pore volume over all cells.
+ double
+ convergenceReduction(const Eigen::Array& B,
+ const Eigen::Array& tempV,
+ const Eigen::Array& R,
+ std::array& R_sum,
+ std::array& maxCoeff,
+ std::array& B_avg,
+ std::vector& maxNormWell,
+ int nc,
+ int nw) const;
+
+ double dpMaxRel() const { return param_.dp_max_rel_; }
+ double dsMax() const { return param_.ds_max_; }
+ double drMaxRel() const { return param_.dr_max_rel_; }
+ double maxResidualAllowed() const { return param_.max_residual_allowed_; }
+
+ };
+} // namespace Opm
+
+#include "BlackoilModel_impl.hpp"
+
+#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED
diff --git a/opm/autodiff/BlackoilModel_impl.hpp b/opm/autodiff/BlackoilModel_impl.hpp
new file mode 100644
index 000000000..1f71307f3
--- /dev/null
+++ b/opm/autodiff/BlackoilModel_impl.hpp
@@ -0,0 +1,2309 @@
+/*
+ Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+ Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
+ Copyright 2014, 2015 Statoil ASA.
+ Copyright 2015 NTNU
+ Copyright 2015 IRIS AS
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
+#define OPM_BLACKOILMODEL_IMPL_HEADER_INCLUDED
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+//#include
+
+// A debugging utility.
+#define OPM_AD_DUMP(foo) \
+ do { \
+ std::cout << "==========================================\n" \
+ << #foo ":\n" \
+ << collapseJacs(foo) << std::endl; \
+ } while (0)
+
+#define OPM_AD_DUMPVAL(foo) \
+ do { \
+ std::cout << "==========================================\n" \
+ << #foo ":\n" \
+ << foo.value() << std::endl; \
+ } while (0)
+
+#define OPM_AD_DISKVAL(foo) \
+ do { \
+ std::ofstream os(#foo); \
+ os.precision(16); \
+ os << foo.value() << std::endl; \
+ } while (0)
+
+
+namespace Opm {
+
+typedef AutoDiffBlock ADB;
+typedef ADB::V V;
+typedef ADB::M M;
+typedef Eigen::Array DataBlock;
+
+
+namespace detail {
+
+
+ std::vector
+ buildAllCells(const int nc)
+ {
+ std::vector all_cells(nc);
+
+ for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
+
+ return all_cells;
+ }
+
+
+
+ template
+ std::vector
+ activePhases(const PU& pu)
+ {
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ std::vector active(maxnp, false);
+
+ for (int p = 0; p < pu.MaxNumPhases; ++p) {
+ active[ p ] = pu.phase_used[ p ] != 0;
+ }
+
+ return active;
+ }
+
+
+
+ template
+ std::vector
+ active2Canonical(const PU& pu)
+ {
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ std::vector act2can(maxnp, -1);
+
+ for (int phase = 0; phase < maxnp; ++phase) {
+ if (pu.phase_used[ phase ]) {
+ act2can[ pu.phase_pos[ phase ] ] = phase;
+ }
+ }
+
+ return act2can;
+ }
+
+
+} // namespace detail
+
+ template
+ void BlackoilModel::ModelParameters::
+ reset()
+ {
+ // default values for the solver parameters
+ dp_max_rel_ = 1.0e9;
+ ds_max_ = 0.2;
+ dr_max_rel_ = 1.0e9;
+ max_residual_allowed_ = 1e7;
+ tolerance_mb_ = 1.0e-5;
+ tolerance_cnv_ = 1.0e-2;
+ tolerance_wells_ = 5.0e-1;
+ }
+
+ template
+ BlackoilModel::ModelParameters::
+ ModelParameters()
+ {
+ // set default values
+ reset();
+ }
+
+ template
+ BlackoilModel::ModelParameters::
+ ModelParameters( const parameter::ParameterGroup& param )
+ {
+ // set default values
+ reset();
+
+ // overload with given parameters
+ dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
+ ds_max_ = param.getDefault("ds_max", ds_max_);
+ dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_);
+ max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_);
+ tolerance_mb_ = param.getDefault("tolerance_mb", tolerance_mb_);
+ tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_);
+ tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
+ }
+
+
+ template
+ BlackoilModel::
+ BlackoilModel(const ModelParameters& param,
+ const Grid& grid ,
+ const BlackoilPropsAdInterface& fluid,
+ const DerivedGeology& geo ,
+ const RockCompressibility* rock_comp_props,
+ const Wells* wells,
+ const NewtonIterationBlackoilInterface& linsolver,
+ const bool has_disgas,
+ const bool has_vapoil,
+ const bool terminal_output)
+ : grid_ (grid)
+ , fluid_ (fluid)
+ , geo_ (geo)
+ , rock_comp_props_(rock_comp_props)
+ , wells_ (wells)
+ , linsolver_ (linsolver)
+ , active_(detail::activePhases(fluid.phaseUsage()))
+ , canph_ (detail::active2Canonical(fluid.phaseUsage()))
+ , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
+ , ops_ (grid)
+ , wops_ (wells_)
+ , has_disgas_(has_disgas)
+ , has_vapoil_(has_vapoil)
+ , param_( param )
+ , use_threshold_pressure_(false)
+ , rq_ (fluid.numPhases())
+ , phaseCondition_(AutoDiffGrid::numCells(grid))
+ , residual_ ( { std::vector(fluid.numPhases(), ADB::null()),
+ ADB::null(),
+ ADB::null() } )
+ , terminal_output_ (terminal_output)
+ {
+#if HAVE_MPI
+ if ( terminal_output_ ) {
+ if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+ {
+ const ParallelISTLInformation& info =
+ boost::any_cast(linsolver_.parallelInformation());
+ // Only rank 0 does print to std::cout if terminal_output is enabled
+ terminal_output_ = (info.communicator().rank()==0);
+ }
+ }
+#endif
+ }
+
+
+
+ template
+ void
+ BlackoilModel::
+ prepareStep(const double dt,
+ ReservoirState& reservoir_state,
+ WellState& /* well_state */)
+ {
+ pvdt_ = geo_.poreVolume() / dt;
+ if (active_[Gas]) {
+ updatePrimalVariableFromState(reservoir_state);
+ }
+ }
+
+
+
+
+ template
+ void
+ BlackoilModel::
+ afterStep(const double /* dt */,
+ ReservoirState& /* reservoir_state */,
+ WellState& /* well_state */)
+ {
+ // Does nothing in this model.
+ }
+
+
+
+
+ template
+ int
+ BlackoilModel::
+ sizeNonLinear() const
+ {
+ return residual_.sizeNonLinear();
+ }
+
+
+
+
+ template
+ int
+ BlackoilModel::
+ linearIterationsLastSolve() const
+ {
+ return linsolver_.iterations();
+ }
+
+
+
+
+ template
+ bool
+ BlackoilModel::
+ terminalOutputEnabled() const
+ {
+ return terminal_output_;
+ }
+
+
+
+
+ template
+ int
+ BlackoilModel::
+ numPhases() const
+ {
+ return fluid_.numPhases();
+ }
+
+
+
+
+ template
+ void
+ BlackoilModel::
+ 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
+ BlackoilModel::ReservoirResidualQuant::ReservoirResidualQuant()
+ : accum(2, ADB::null())
+ , mflux( ADB::null())
+ , b ( ADB::null())
+ , head ( ADB::null())
+ , mob ( ADB::null())
+ {
+ }
+
+
+
+
+
+ template
+ BlackoilModel::SolutionState::SolutionState(const int np)
+ : pressure ( ADB::null())
+ , temperature( ADB::null())
+ , saturation(np, ADB::null())
+ , rs ( ADB::null())
+ , rv ( ADB::null())
+ , qs ( ADB::null())
+ , bhp ( ADB::null())
+ , canonical_phase_pressures(3, ADB::null())
+ {
+ }
+
+
+
+
+
+ template
+ BlackoilModel::
+ WellOps::WellOps(const Wells* wells)
+ : w2p(),
+ p2w()
+ {
+ if( wells )
+ {
+ w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
+ p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
+
+ const int nw = wells->number_of_wells;
+ const int* const wpos = wells->well_connpos;
+
+ typedef Eigen::Triplet Tri;
+
+ std::vector scatter, gather;
+ scatter.reserve(wpos[nw]);
+ gather .reserve(wpos[nw]);
+
+ for (int w = 0, i = 0; w < nw; ++w) {
+ for (; i < wpos[ w + 1 ]; ++i) {
+ scatter.push_back(Tri(i, w, 1.0));
+ gather .push_back(Tri(w, i, 1.0));
+ }
+ }
+
+ w2p.setFromTriplets(scatter.begin(), scatter.end());
+ p2w.setFromTriplets(gather .begin(), gather .end());
+ }
+ }
+
+
+
+
+
+ template
+ typename BlackoilModel::SolutionState
+ BlackoilModel::constantState(const BlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw) const
+ {
+ auto state = variableState(x, xw);
+ makeConstantState(state);
+ return state;
+ }
+
+
+
+
+
+ template
+ void
+ BlackoilModel::makeConstantState(SolutionState& state) const
+ {
+ // HACK: throw away the derivatives. this may not be the most
+ // performant way to do things, but it will make the state
+ // automatically consistent with variableState() (and doing
+ // things automatically is all the rage in this module ;)
+ state.pressure = ADB::constant(state.pressure.value());
+ state.temperature = ADB::constant(state.temperature.value());
+ state.rs = ADB::constant(state.rs.value());
+ state.rv = ADB::constant(state.rv.value());
+ const int num_phases = state.saturation.size();
+ for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
+ state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
+ }
+ state.qs = ADB::constant(state.qs.value());
+ state.bhp = ADB::constant(state.bhp.value());
+ assert(state.canonical_phase_pressures.size() == static_cast(Opm::BlackoilPhases::MaxNumPhases));
+ for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) {
+ ADB& pp = state.canonical_phase_pressures[canphase];
+ pp = ADB::constant(pp.value());
+ }
+ }
+
+
+
+
+
+ template
+ typename BlackoilModel::SolutionState
+ BlackoilModel::variableState(const BlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw) const
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+ const int np = x.numPhases();
+
+ std::vector vars0;
+ // p, Sw and Rs, Rv or Sg is used as primary depending on solution conditions
+ vars0.reserve(np + 1);
+ // Initial pressure.
+ assert (not x.pressure().empty());
+ const V p = Eigen::Map(& x.pressure()[0], nc, 1);
+ vars0.push_back(p);
+
+ // Initial saturation.
+ assert (not x.saturation().empty());
+ const DataBlock s = Eigen::Map(& x.saturation()[0], nc, np);
+ const Opm::PhaseUsage pu = fluid_.phaseUsage();
+ // We do not handle a Water/Gas situation correctly, guard against it.
+ assert (active_[ Oil]);
+ if (active_[ Water ]) {
+ const V sw = s.col(pu.phase_pos[ Water ]);
+ vars0.push_back(sw);
+ }
+
+ // store cell status in vectors
+ V isRs = V::Zero(nc,1);
+ V isRv = V::Zero(nc,1);
+ V isSg = V::Zero(nc,1);
+
+ if (active_[ Gas ]){
+ for (int c = 0; c < nc ; c++ ) {
+ switch (primalVariable_[c]) {
+ case PrimalVariables::RS:
+ isRs[c] = 1;
+ break;
+
+ case PrimalVariables::RV:
+ isRv[c] = 1;
+ break;
+
+ default:
+ isSg[c] = 1;
+ break;
+ }
+ }
+
+
+ // define new primary variable xvar depending on solution condition
+ V xvar(nc);
+ const V sg = s.col(pu.phase_pos[ Gas ]);
+ const V rs = Eigen::Map(& x.gasoilratio()[0], x.gasoilratio().size());
+ const V rv = Eigen::Map(& x.rv()[0], x.rv().size());
+ xvar = isRs*rs + isRv*rv + isSg*sg;
+ vars0.push_back(xvar);
+ }
+
+
+ // Initial well rates.
+ if( wellsActive() )
+ {
+ // Need to reshuffle well rates, from phase running fastest
+ // to wells running fastest.
+ const int nw = wells().number_of_wells;
+
+ // The transpose() below switches the ordering.
+ const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose();
+ const V qs = Eigen::Map(wrates.data(), nw*np);
+ vars0.push_back(qs);
+
+ // Initial well bottom-hole pressure.
+ assert (not xw.bhp().empty());
+ const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size());
+ vars0.push_back(bhp);
+ }
+ else
+ {
+ // push null states for qs and bhp
+ vars0.push_back(V());
+ vars0.push_back(V());
+ }
+
+ std::vector vars = ADB::variables(vars0);
+
+ SolutionState state(np);
+
+ // Pressure.
+ int nextvar = 0;
+ state.pressure = std::move(vars[ nextvar++ ]);
+
+ // temperature
+ const V temp = Eigen::Map(& x.temperature()[0], x.temperature().size());
+ state.temperature = ADB::constant(temp);
+
+ // Saturations
+ {
+
+ ADB so = ADB::constant(V::Ones(nc, 1));
+
+ if (active_[ Water ]) {
+ state.saturation[pu.phase_pos[ Water ]] = std::move(vars[ nextvar++ ]);
+ const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
+ so -= sw;
+ }
+
+ if (active_[ Gas ]) {
+ // Define Sg Rs and Rv in terms of xvar.
+ // Xvar is only defined if gas phase is active
+ const ADB& xvar = vars[ nextvar++ ];
+ ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
+ sg = isSg*xvar + isRv* so;
+ so -= sg;
+
+ if (active_[ Oil ]) {
+ // RS and RV is only defined if both oil and gas phase are active.
+ const ADB& sw = (active_[ Water ]
+ ? state.saturation[ pu.phase_pos[ Water ] ]
+ : ADB::constant(V::Zero(nc, 1)));
+ state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
+ const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
+ if (has_disgas_) {
+ state.rs = (1-isRs) * rsSat + isRs*xvar;
+ } else {
+ state.rs = rsSat;
+ }
+ const ADB rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
+ if (has_vapoil_) {
+ state.rv = (1-isRv) * rvSat + isRv*xvar;
+ } else {
+ state.rv = rvSat;
+ }
+ }
+ }
+
+ if (active_[ Oil ]) {
+ // Note that so is never a primary variable.
+ state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
+ }
+ }
+
+ // Qs.
+ state.qs = std::move(vars[ nextvar++ ]);
+
+ // Bhp.
+ state.bhp = std::move(vars[ nextvar++ ]);
+
+ assert(nextvar == int(vars.size()));
+
+ return state;
+ }
+
+
+
+
+
+ template
+ void
+ BlackoilModel::computeAccum(const SolutionState& state,
+ const int aix )
+ {
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+ const ADB& press = state.pressure;
+ const ADB& temp = state.temperature;
+ const std::vector& sat = state.saturation;
+ const ADB& rs = state.rs;
+ const ADB& rv = state.rv;
+
+ const std::vector cond = phaseCondition();
+
+ 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, state.canonical_phase_pressures[phase], temp, rs, rv, cond, cells_);
+ rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
+ // OPM_AD_DUMP(rq_[pos].b);
+ // OPM_AD_DUMP(rq_[pos].accum[aix]);
+ }
+ }
+
+ if (active_[ Oil ] && active_[ Gas ]) {
+ // Account for gas dissolved in oil and vaporized oil
+ const int po = pu.phase_pos[ Oil ];
+ const int pg = pu.phase_pos[ Gas ];
+
+ // Temporary copy to avoid contribution of dissolved gas in the vaporized oil
+ // when both dissolved gas and vaporized oil are present.
+ const ADB accum_gas_copy =rq_[pg].accum[aix];
+
+ rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix];
+ rq_[po].accum[aix] += state.rv * accum_gas_copy;
+ // OPM_AD_DUMP(rq_[pg].accum[aix]);
+ }
+ }
+
+
+
+
+
+ template
+ void BlackoilModel::computeWellConnectionPressures(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw)
+ {
+ if( ! wellsActive() ) return ;
+
+ using namespace Opm::AutoDiffGrid;
+ // 1. Compute properties required by computeConnectionPressureDelta().
+ // Note that some of the complexity of this part is due to the function
+ // taking std::vector arguments, and not Eigen objects.
+ const int nperf = wells().well_connpos[wells().number_of_wells];
+ const int nw = wells().number_of_wells;
+ const std::vector well_cells(wells().well_cells, wells().well_cells + nperf);
+
+ // Compute the average pressure in each well block
+ const V perf_press = Eigen::Map(xw.perfPress().data(), nperf);
+ V avg_press = perf_press*0;
+ for (int w = 0; w < nw; ++w) {
+ for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
+ const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
+ const double p_avg = (perf_press[perf] + p_above)/2;
+ avg_press[perf] = p_avg;
+ }
+ }
+
+ // Use cell values for the temperature as the wells don't knows its temperature yet.
+ const ADB perf_temp = subset(state.temperature, well_cells);
+
+ // Compute b, rsmax, rvmax values for perforations.
+ // Evaluate the properties using average well block pressures
+ // and cell values for rs, rv, phase condition and temperature.
+ const ADB avg_press_ad = ADB::constant(avg_press);
+ std::vector perf_cond(nperf);
+ const std::vector& pc = phaseCondition();
+ for (int perf = 0; perf < nperf; ++perf) {
+ perf_cond[perf] = pc[well_cells[perf]];
+ }
+ const PhaseUsage& pu = fluid_.phaseUsage();
+ DataBlock b(nperf, pu.num_phases);
+ std::vector rsmax_perf(nperf, 0.0);
+ std::vector rvmax_perf(nperf, 0.0);
+ if (pu.phase_used[BlackoilPhases::Aqua]) {
+ const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
+ b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
+ }
+ assert(active_[Oil]);
+ const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
+ if (pu.phase_used[BlackoilPhases::Liquid]) {
+ const ADB perf_rs = subset(state.rs, well_cells);
+ const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
+ b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
+ const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
+ rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
+ }
+ if (pu.phase_used[BlackoilPhases::Vapour]) {
+ const ADB perf_rv = subset(state.rv, well_cells);
+ const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
+ b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
+ const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
+ rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
+ }
+ // b is row major, so can just copy data.
+ std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases);
+ // Extract well connection depths.
+ const V depth = cellCentroidsZToEigen(grid_);
+ const V pdepth = subset(depth, well_cells);
+ std::vector perf_depth(pdepth.data(), pdepth.data() + nperf);
+ // Surface density.
+ std::vector surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
+ // Gravity
+ double grav = 0.0;
+ const double* g = geo_.gravity();
+ const int dim = dimensions(grid_);
+ if (g) {
+ // Guard against gravity in anything but last dimension.
+ for (int dd = 0; dd < dim - 1; ++dd) {
+ assert(g[dd] == 0.0);
+ }
+ grav = g[dim - 1];
+ }
+
+ // 2. Compute pressure deltas, and store the results.
+ std::vector cdp = WellDensitySegmented
+ ::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(),
+ b_perf, rsmax_perf, rvmax_perf, perf_depth,
+ surf_dens, grav);
+ well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf);
+ }
+
+
+
+
+
+ template
+ void
+ BlackoilModel::
+ assemble(const BlackoilState& reservoir_state,
+ WellStateFullyImplicitBlackoil& well_state,
+ const bool initial_assembly)
+ {
+ using namespace Opm::AutoDiffGrid;
+
+ // Possibly switch well controls and updating well state to
+ // get reasonable initial conditions for the wells
+ updateWellControls(well_state);
+
+ // Create the primary variables.
+ SolutionState state = variableState(reservoir_state, well_state);
+
+ if (initial_assembly) {
+ // Create the (constant, derivativeless) initial state.
+ SolutionState state0 = state;
+ makeConstantState(state0);
+ // Compute initial accumulation contributions
+ // and well connection pressures.
+ computeAccum(state0, 0);
+ computeWellConnectionPressures(state0, well_state);
+ }
+
+ // OPM_AD_DISKVAL(state.pressure);
+ // OPM_AD_DISKVAL(state.saturation[0]);
+ // OPM_AD_DISKVAL(state.saturation[1]);
+ // OPM_AD_DISKVAL(state.saturation[2]);
+ // OPM_AD_DISKVAL(state.rs);
+ // OPM_AD_DISKVAL(state.rv);
+ // OPM_AD_DISKVAL(state.qs);
+ // OPM_AD_DISKVAL(state.bhp);
+
+ // -------- Mass balance equations --------
+
+ // Compute b_p and the accumulation term b_p*s_p for each phase,
+ // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
+ // These quantities are stored in rq_[phase].accum[1].
+ // The corresponding accumulation terms from the start of
+ // the timestep (b^0_p*s^0_p etc.) were already computed
+ // on the initial call to assemble() and stored in rq_[phase].accum[0].
+ computeAccum(state, 1);
+
+ // Set up the common parts of the mass balance equations
+ // for each active phase.
+ const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
+ const std::vector kr = computeRelPerm(state);
+ for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
+ computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
+
+ residual_.material_balance_eq[ phaseIdx ] =
+ pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ + ops_.div*rq_[phaseIdx].mflux;
+ }
+
+ // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
+
+ // Add the extra (flux) terms to the mass balance equations
+ // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
+ // The extra terms in the accumulation part of the equation are already handled.
+ if (active_[ Oil ] && active_[ Gas ]) {
+ const int po = fluid_.phaseUsage().phase_pos[ Oil ];
+ const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
+
+ const UpwindSelector upwindOil(grid_, ops_,
+ rq_[po].head.value());
+ const ADB rs_face = upwindOil.select(state.rs);
+
+ const UpwindSelector upwindGas(grid_, ops_,
+ rq_[pg].head.value());
+ const ADB rv_face = upwindGas.select(state.rv);
+
+ residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
+ residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
+
+ // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
+
+ }
+
+ // -------- Well equations ----------
+
+ // Add contribution from wells and set up the well equations.
+ V aliveWells;
+ addWellEq(state, well_state, aliveWells);
+ addWellControlEq(state, well_state, aliveWells);
+ }
+
+
+
+
+
+ template
+ void BlackoilModel::addWellEq(const SolutionState& state,
+ WellStateFullyImplicitBlackoil& xw,
+ V& aliveWells)
+ {
+ if( ! wellsActive() ) return ;
+
+ const int nc = Opm::AutoDiffGrid::numCells(grid_);
+ const int np = wells().number_of_phases;
+ const int nw = wells().number_of_wells;
+ const int nperf = wells().well_connpos[nw];
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ V Tw = Eigen::Map(wells().WI, nperf);
+ const std::vector well_cells(wells().well_cells, wells().well_cells + nperf);
+
+ // pressure diffs computed already (once per step, not changing per iteration)
+ const V& cdp = well_perforation_pressure_diffs_;
+
+ // Extract needed quantities for the perforation cells
+ const ADB& p_perfcells = subset(state.pressure, well_cells);
+ const ADB& rv_perfcells = subset(state.rv,well_cells);
+ const ADB& rs_perfcells = subset(state.rs,well_cells);
+ std::vector mob_perfcells(np, ADB::null());
+ std::vector b_perfcells(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ mob_perfcells[phase] = subset(rq_[phase].mob,well_cells);
+ b_perfcells[phase] = subset(rq_[phase].b,well_cells);
+ }
+
+ // Perforation pressure
+ const ADB perfpressure = (wops_.w2p * state.bhp) + cdp;
+ std::vector perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf);
+ xw.perfPress() = perfpressure_d;
+
+ // Pressure drawdown (also used to determine direction of flow)
+ const ADB drawdown = p_perfcells - perfpressure;
+
+ // Compute vectors with zero and ones that
+ // selects the wanted quantities.
+
+ // selects injection perforations
+ V selectInjectingPerforations = V::Zero(nperf);
+ // selects producing perforations
+ V selectProducingPerforations = V::Zero(nperf);
+ for (int c = 0; c < nperf; ++c){
+ if (drawdown.value()[c] < 0)
+ selectInjectingPerforations[c] = 1;
+ else
+ selectProducingPerforations[c] = 1;
+ }
+
+ // HANDLE FLOW INTO WELLBORE
+
+ // compute phase volumetric rates at standard conditions
+ std::vector cq_ps(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
+ cq_ps[phase] = b_perfcells[phase] * cq_p;
+ }
+ if (active_[Oil] && active_[Gas]) {
+ const int oilpos = pu.phase_pos[Oil];
+ const int gaspos = pu.phase_pos[Gas];
+ const ADB cq_psOil = cq_ps[oilpos];
+ const ADB cq_psGas = cq_ps[gaspos];
+ cq_ps[gaspos] += rs_perfcells * cq_psOil;
+ cq_ps[oilpos] += rv_perfcells * cq_psGas;
+ }
+
+ // HANDLE FLOW OUT FROM WELLBORE
+
+ // Using total mobilities
+ ADB total_mob = mob_perfcells[0];
+ for (int phase = 1; phase < np; ++phase) {
+ total_mob += mob_perfcells[phase];
+ }
+ // injection perforations total volume rates
+ const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
+
+ // compute wellbore mixture for injecting perforations
+ // The wellbore mixture depends on the inflow from the reservoar
+ // and the well injection rates.
+
+ // compute avg. and total wellbore phase volumetric rates at standard conds
+ const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np);
+ std::vector wbq(np, ADB::null());
+ ADB wbqt = ADB::constant(V::Zero(nw));
+ for (int phase = 0; phase < np; ++phase) {
+ const ADB& q_ps = wops_.p2w * cq_ps[phase];
+ const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
+ Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero);
+ const int pos = pu.phase_pos[phase];
+ wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps;
+ wbqt += wbq[phase];
+ }
+
+ // compute wellbore mixture at standard conditions.
+ Selector notDeadWells_selector(wbqt.value(), Selector::Zero);
+ std::vector cmix_s(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ const int pos = pu.phase_pos[phase];
+ cmix_s[phase] = wops_.w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
+ }
+
+ // compute volume ratio between connection at standard conditions
+ ADB volumeRatio = ADB::constant(V::Zero(nperf));
+ const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
+ for (int phase = 0; phase < np; ++phase) {
+ ADB tmp = cmix_s[phase];
+
+ if (phase == Oil && active_[Gas]) {
+ const int gaspos = pu.phase_pos[Gas];
+ tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
+ }
+ if (phase == Gas && active_[Oil]) {
+ const int oilpos = pu.phase_pos[Oil];
+ tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
+ }
+ volumeRatio += tmp / b_perfcells[phase];
+ }
+
+ // injecting connections total volumerates at standard conditions
+ ADB cqt_is = cqt_i/volumeRatio;
+
+ // connection phase volumerates at standard conditions
+ std::vector cq_s(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
+ }
+
+ // Add well contributions to mass balance equations
+ for (int phase = 0; phase < np; ++phase) {
+ residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
+ }
+
+ // WELL EQUATIONS
+ ADB qs = state.qs;
+ for (int phase = 0; phase < np; ++phase) {
+ qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
+
+ }
+
+ // check for dead wells (used in the well controll equations)
+ aliveWells = V::Constant(nw, 1.0);
+ for (int w = 0; w < nw; ++w) {
+ if (wbqt.value()[w] == 0) {
+ aliveWells[w] = 0.0;
+ }
+ }
+
+ // Update the perforation phase rates (used to calculate the pressure drop in the wellbore)
+ V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
+ for (int phase = 1; phase < np; ++phase) {
+ cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
+ }
+
+ std::vector cq_d(cq.data(), cq.data() + nperf*np);
+ xw.perfPhaseRates() = cq_d;
+
+ residual_.well_flux_eq = qs;
+ }
+
+
+
+
+
+ namespace detail
+ {
+ double rateToCompare(const std::vector& well_phase_flow_rate,
+ const int well,
+ const int num_phases,
+ const double* distr)
+ {
+ double rate = 0.0;
+ for (int phase = 0; phase < num_phases; ++phase) {
+ // Important: well_phase_flow_rate is ordered with all phase rates for first
+ // well first, then all phase rates for second well etc.
+ rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
+ }
+ return rate;
+ }
+
+ bool constraintBroken(const std::vector& bhp,
+ const std::vector& well_phase_flow_rate,
+ const int well,
+ const int num_phases,
+ const WellType& well_type,
+ const WellControls* wc,
+ const int ctrl_index)
+ {
+ const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
+ const double target = well_controls_iget_target(wc, ctrl_index);
+ const double* distr = well_controls_iget_distr(wc, ctrl_index);
+
+ bool broken = false;
+
+ switch (well_type) {
+ case INJECTOR:
+ {
+ switch (ctrl_type) {
+ case BHP:
+ broken = bhp[well] > target;
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ broken = rateToCompare(well_phase_flow_rate,
+ well, num_phases, distr) > target;
+ break;
+ }
+ }
+ break;
+
+ case PRODUCER:
+ {
+ switch (ctrl_type) {
+ case BHP:
+ broken = bhp[well] < target;
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ // Note that the rates compared below are negative,
+ // so breaking the constraints means: too high flow rate
+ // (as for injection).
+ broken = rateToCompare(well_phase_flow_rate,
+ well, num_phases, distr) < target;
+ break;
+ }
+ }
+ break;
+
+ default:
+ OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
+ }
+
+ return broken;
+ }
+ } // namespace detail
+
+
+
+
+ template
+ void BlackoilModel::updateWellControls(WellStateFullyImplicitBlackoil& xw) const
+ {
+ if( ! wellsActive() ) return ;
+
+ std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
+ // Find, for each well, if any constraints are broken. If so,
+ // switch control to first broken constraint.
+ const int np = wells().number_of_phases;
+ const int nw = wells().number_of_wells;
+ for (int w = 0; w < nw; ++w) {
+ const WellControls* wc = wells().ctrls[w];
+ // The current control in the well state overrides
+ // the current control set in the Wells struct, which
+ // is instead treated as a default.
+ int current = xw.currentControls()[w];
+ // Loop over all controls except the current one, and also
+ // skip any RESERVOIR_RATE controls, since we cannot
+ // handle those.
+ const int nwc = well_controls_get_num(wc);
+ int ctrl_index = 0;
+ for (; ctrl_index < nwc; ++ctrl_index) {
+ if (ctrl_index == current) {
+ // This is the currently used control, so it is
+ // used as an equation. So this is not used as an
+ // inequality constraint, and therefore skipped.
+ continue;
+ }
+ if (detail::constraintBroken(xw.bhp(), xw.wellRates(), w, np, wells().type[w], wc, ctrl_index)) {
+ // ctrl_index will be the index of the broken constraint after the loop.
+ break;
+ }
+ }
+ if (ctrl_index != nwc) {
+ // Constraint number ctrl_index was broken, switch to it.
+ if (terminal_output_)
+ {
+ std::cout << "Switching control mode for well " << wells().name[w]
+ << " from " << modestring[well_controls_iget_type(wc, current)]
+ << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
+ }
+ xw.currentControls()[w] = ctrl_index;
+ current = xw.currentControls()[w];
+ }
+
+ // Updating well state and primary variables.
+ // Target values are used as initial conditions for BHP and SURFACE_RATE
+ const double target = well_controls_iget_target(wc, current);
+ const double* distr = well_controls_iget_distr(wc, current);
+ switch (well_controls_iget_type(wc, current)) {
+ case BHP:
+ xw.bhp()[w] = target;
+ break;
+
+ case RESERVOIR_RATE:
+ // No direct change to any observable quantity at
+ // surface condition. In this case, use existing
+ // flow rates as initial conditions as reservoir
+ // rate acts only in aggregate.
+ break;
+
+ case SURFACE_RATE:
+ for (int phase = 0; phase < np; ++phase) {
+ if (distr[phase] > 0.0) {
+ xw.wellRates()[np*w + phase] = target * distr[phase];
+ }
+ }
+ break;
+ }
+
+ }
+ }
+
+
+
+
+ template
+ void BlackoilModel::addWellControlEq(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw,
+ const V& aliveWells)
+ {
+ if( ! wellsActive() ) return;
+
+ const int np = wells().number_of_phases;
+ const int nw = wells().number_of_wells;
+
+ V bhp_targets = V::Zero(nw);
+ V rate_targets = V::Zero(nw);
+ M rate_distr(nw, np*nw);
+ for (int w = 0; w < nw; ++w) {
+ const WellControls* wc = wells().ctrls[w];
+ // The current control in the well state overrides
+ // the current control set in the Wells struct, which
+ // is instead treated as a default.
+ const int current = xw.currentControls()[w];
+
+ switch (well_controls_iget_type(wc, current)) {
+ case BHP:
+ {
+ bhp_targets (w) = well_controls_iget_target(wc, current);
+ rate_targets(w) = -1e100;
+ }
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ {
+ // RESERVOIR and SURFACE rates look the same, from a
+ // high-level point of view, in the system of
+ // simultaneous linear equations.
+
+ const double* const distr =
+ well_controls_iget_distr(wc, current);
+
+ for (int p = 0; p < np; ++p) {
+ rate_distr.insert(w, p*nw + w) = distr[p];
+ }
+
+ bhp_targets (w) = -1.0e100;
+ rate_targets(w) = well_controls_iget_target(wc, current);
+ }
+ break;
+ }
+ }
+ const ADB bhp_residual = state.bhp - bhp_targets;
+ const ADB rate_residual = rate_distr * state.qs - rate_targets;
+ // Choose bhp residual for positive bhp targets.
+ Selector bhp_selector(bhp_targets);
+ residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
+ // For wells that are dead (not flowing), and therefore not communicating
+ // with the reservoir, we set the equation to be equal to the well's total
+ // flow. This will be a solution only if the target rate is also zero.
+ M rate_summer(nw, np*nw);
+ for (int w = 0; w < nw; ++w) {
+ for (int phase = 0; phase < np; ++phase) {
+ rate_summer.insert(w, phase*nw + w) = 1.0;
+ }
+ }
+ Selector alive_selector(aliveWells, Selector::NotEqualZero);
+ residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs);
+ // OPM_AD_DUMP(residual_.well_eq);
+ }
+
+
+
+
+
+ template
+ V BlackoilModel::solveJacobianSystem() const
+ {
+ return linsolver_.computeNewtonIncrement(residual_);
+ }
+
+
+
+
+
+ namespace detail
+ {
+ /// \brief Compute the L-infinity norm of a vector
+ /// \warn This function is not suitable to compute on the well equations.
+ /// \param a The container to compute the infinity norm on.
+ /// It has to have one entry for each cell.
+ /// \param info In a parallel this holds the information about the data distribution.
+ double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() )
+ {
+#if HAVE_MPI
+ if ( pinfo.type() == typeid(ParallelISTLInformation) )
+ {
+ const ParallelISTLInformation& real_info =
+ boost::any_cast(pinfo);
+ double result=0;
+ real_info.computeReduction(a.value(), Reduction::makeGlobalMaxFunctor(), result);
+ return result;
+ }
+ else
+#endif
+ {
+ if( a.value().size() > 0 ) {
+ return a.value().matrix().lpNorm ();
+ }
+ else { // this situation can occur when no wells are present
+ return 0.0;
+ }
+ }
+ }
+
+ /// \brief Compute the L-infinity norm of a vector representing a well equation.
+ /// \param a The container to compute the infinity norm on.
+ /// \param info In a parallel this holds the information about the data distribution.
+ double infinityNormWell( const ADB& a, const boost::any& pinfo )
+ {
+ double result=0;
+ if( a.value().size() > 0 ) {
+ result = a.value().matrix().lpNorm ();
+ }
+#if HAVE_MPI
+ if ( pinfo.type() == typeid(ParallelISTLInformation) )
+ {
+ const ParallelISTLInformation& real_info =
+ boost::any_cast(pinfo);
+ result = real_info.communicator().max(result);
+ }
+#endif
+ return result;
+ }
+
+ } // namespace detail
+
+
+
+
+
+ template
+ void BlackoilModel::updateState(const V& dx,
+ BlackoilState& reservoir_state,
+ WellStateFullyImplicitBlackoil& well_state)
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int np = fluid_.numPhases();
+ const int nc = numCells(grid_);
+ const int nw = wellsActive() ? wells().number_of_wells : 0;
+ const V null;
+ assert(null.size() == 0);
+ const V zero = V::Zero(nc);
+
+ // store cell status in vectors
+ V isRs = V::Zero(nc,1);
+ V isRv = V::Zero(nc,1);
+ V isSg = V::Zero(nc,1);
+ if (active_[Gas]) {
+ for (int c = 0; c < nc; ++c) {
+ switch (primalVariable_[c]) {
+ case PrimalVariables::RS:
+ isRs[c] = 1;
+ break;
+
+ case PrimalVariables::RV:
+ isRv[c] = 1;
+ break;
+
+ default:
+ isSg[c] = 1;
+ break;
+ }
+ }
+ }
+
+ // Extract parts of dx corresponding to each part.
+ const V dp = subset(dx, Span(nc));
+ int varstart = nc;
+ const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
+ varstart += dsw.size();
+
+ const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
+ varstart += dxvar.size();
+
+ const V dqs = subset(dx, Span(np*nw, 1, varstart));
+ varstart += dqs.size();
+ const V dbhp = subset(dx, Span(nw, 1, varstart));
+ varstart += dbhp.size();
+ assert(varstart == dx.size());
+
+ // Pressure update.
+ const double dpmaxrel = dpMaxRel();
+ const V p_old = Eigen::Map(&reservoir_state.pressure()[0], nc, 1);
+ const V absdpmax = dpmaxrel*p_old.abs();
+ const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
+ const V p = (p_old - dp_limited).max(zero);
+ std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
+
+
+ // Saturation updates.
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const DataBlock s_old = Eigen::Map(& reservoir_state.saturation()[0], nc, np);
+ const double dsmax = dsMax();
+
+ V so;
+ V sw;
+ V sg;
+
+ {
+ V maxVal = zero;
+ V dso = zero;
+ if (active_[Water]){
+ maxVal = dsw.abs().max(maxVal);
+ dso = dso - dsw;
+ }
+
+ V dsg;
+ if (active_[Gas]){
+ dsg = isSg * dxvar - isRv * dsw;
+ maxVal = dsg.abs().max(maxVal);
+ dso = dso - dsg;
+ }
+
+ maxVal = dso.abs().max(maxVal);
+
+ V step = dsmax/maxVal;
+ step = step.min(1.);
+
+ if (active_[Water]) {
+ const int pos = pu.phase_pos[ Water ];
+ const V sw_old = s_old.col(pos);
+ sw = sw_old - step * dsw;
+ }
+
+ if (active_[Gas]) {
+ const int pos = pu.phase_pos[ Gas ];
+ const V sg_old = s_old.col(pos);
+ sg = sg_old - step * dsg;
+ }
+
+ const int pos = pu.phase_pos[ Oil ];
+ const V so_old = s_old.col(pos);
+ so = so_old - step * dso;
+ }
+
+ // Appleyard chop process.
+ auto ixg = sg < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixg[c]) {
+ sw[c] = sw[c] / (1-sg[c]);
+ so[c] = so[c] / (1-sg[c]);
+ sg[c] = 0;
+ }
+ }
+
+
+ auto ixo = so < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixo[c]) {
+ sw[c] = sw[c] / (1-so[c]);
+ sg[c] = sg[c] / (1-so[c]);
+ so[c] = 0;
+ }
+ }
+
+ auto ixw = sw < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixw[c]) {
+ so[c] = so[c] / (1-sw[c]);
+ sg[c] = sg[c] / (1-sw[c]);
+ sw[c] = 0;
+ }
+ }
+
+ const V sumSat = sw + so + sg;
+ sw = sw / sumSat;
+ so = so / sumSat;
+ sg = sg / sumSat;
+
+ // Update the reservoir_state
+ for (int c = 0; c < nc; ++c) {
+ reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
+ }
+
+ for (int c = 0; c < nc; ++c) {
+ reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
+ }
+
+ if (active_[ Oil ]) {
+ const int pos = pu.phase_pos[ Oil ];
+ for (int c = 0; c < nc; ++c) {
+ reservoir_state.saturation()[c*np + pos] = so[c];
+ }
+ }
+
+ // Update rs and rv
+ const double drmaxrel = drMaxRel();
+ V rs;
+ if (has_disgas_) {
+ const V rs_old = Eigen::Map(&reservoir_state.gasoilratio()[0], nc);
+ const V drs = isRs * dxvar;
+ const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
+ rs = rs_old - drs_limited;
+ }
+ V rv;
+ if (has_vapoil_) {
+ const V rv_old = Eigen::Map(&reservoir_state.rv()[0], nc);
+ const V drv = isRv * dxvar;
+ const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
+ rv = rv_old - drv_limited;
+ }
+
+
+ // Sg is used as primal variable for water only cells.
+ const double epsilon = std::sqrt(std::numeric_limits::epsilon());
+ auto watOnly = sw > (1 - epsilon);
+
+ // phase translation sg <-> rs
+ std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
+
+ if (has_disgas_) {
+ const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
+ const V rsSat = fluidRsSat(p, so, cells_);
+ // The obvious case
+ auto hasGas = (sg > 0 && isRs == 0);
+
+ // Set oil saturated if previous rs is sufficiently large
+ const V rs_old = Eigen::Map(&reservoir_state.gasoilratio()[0], nc);
+ auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
+ auto useSg = watOnly || hasGas || gasVaporized;
+ for (int c = 0; c < nc; ++c) {
+ if (useSg[c]) {
+ rs[c] = rsSat[c];
+ } else {
+ primalVariable_[c] = PrimalVariables::RS;
+ }
+ }
+
+ }
+
+ // phase transitions so <-> rv
+ if (has_vapoil_) {
+
+ // The gas pressure is needed for the rvSat calculations
+ const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
+ const V gaspress = computeGasPressure(p, sw, so, sg);
+ const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
+ const V rvSat = fluidRvSat(gaspress, so, cells_);
+
+ // The obvious case
+ auto hasOil = (so > 0 && isRv == 0);
+
+ // Set oil saturated if previous rv is sufficiently large
+ const V rv_old = Eigen::Map(&reservoir_state.rv()[0], nc);
+ auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) );
+ auto useSg = watOnly || hasOil || oilCondensed;
+ for (int c = 0; c < nc; ++c) {
+ if (useSg[c]) {
+ rv[c] = rvSat[c];
+ } else {
+ primalVariable_[c] = PrimalVariables::RV;
+ }
+ }
+
+ }
+
+ // Update the reservoir_state
+ if (has_disgas_) {
+ std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
+ }
+
+ if (has_vapoil_) {
+ std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
+ }
+
+ if( wellsActive() )
+ {
+ // Qs update.
+ // Since we need to update the wellrates, that are ordered by wells,
+ // from dqs which are ordered by phase, the simplest is to compute
+ // dwr, which is the data from dqs but ordered by wells.
+ const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose();
+ const V dwr = Eigen::Map(wwr.data(), nw*np);
+ const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np);
+ const V wr = wr_old - dwr;
+ std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
+
+ // Bhp update.
+ const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1);
+ const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
+ const V bhp = bhp_old - dbhp_limited;
+ std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
+ }
+
+ // Update phase conditions used for property calculations.
+ updatePhaseCondFromPrimalVariable();
+ }
+
+
+
+
+
+ template
+ std::vector
+ BlackoilModel::computeRelPerm(const SolutionState& state) const
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+
+ const ADB zero = ADB::constant(V::Zero(nc));
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const ADB& sw = (active_[ Water ]
+ ? state.saturation[ pu.phase_pos[ Water ] ]
+ : zero);
+
+ const ADB& so = (active_[ Oil ]
+ ? state.saturation[ pu.phase_pos[ Oil ] ]
+ : zero);
+
+ const ADB& sg = (active_[ Gas ]
+ ? state.saturation[ pu.phase_pos[ Gas ] ]
+ : zero);
+
+ return fluid_.relperm(sw, so, sg, cells_);
+ }
+
+
+
+
+
+ template
+ std::vector
+ BlackoilModel::computePressures(const SolutionState& state) const
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+
+ const ADB null = ADB::constant(V::Zero(nc));
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const ADB& sw = (active_[ Water ]
+ ? state.saturation[ pu.phase_pos[ Water ] ]
+ : null);
+
+ const ADB& so = (active_[ Oil ]
+ ? state.saturation[ pu.phase_pos[ Oil ] ]
+ : null);
+
+ const ADB& sg = (active_[ Gas ]
+ ? state.saturation[ pu.phase_pos[ Gas ] ]
+ : null);
+ return computePressures(state.pressure, sw, so, sg);
+
+
+ }
+
+
+
+
+
+ template
+ std::vector
+ BlackoilModel::
+ computePressures(const ADB& po,
+ const ADB& sw,
+ const ADB& so,
+ const ADB& sg) const
+ {
+ // convert the pressure offsets to the capillary pressures
+ std::vector pressure = fluid_.capPress(sw, so, sg, cells_);
+ for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
+ // The reference pressure is always the liquid phase (oil) pressure.
+ if (phaseIdx == BlackoilPhases::Liquid)
+ continue;
+ pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
+ }
+
+ // Since pcow = po - pw, but pcog = pg - po,
+ // we have
+ // pw = po - pcow
+ // pg = po + pcgo
+ // This is an unfortunate inconsistency, but a convention we must handle.
+ for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
+ if (phaseIdx == BlackoilPhases::Aqua) {
+ pressure[phaseIdx] = po - pressure[phaseIdx];
+ } else {
+ pressure[phaseIdx] += po;
+ }
+ }
+
+ return pressure;
+ }
+
+
+
+
+
+ template
+ V
+ BlackoilModel::computeGasPressure(const V& po,
+ const V& sw,
+ const V& so,
+ const V& sg) const
+ {
+ assert (active_[Gas]);
+ std::vector cp = fluid_.capPress(ADB::constant(sw),
+ ADB::constant(so),
+ ADB::constant(sg),
+ cells_);
+ return cp[Gas].value() + po;
+ }
+
+
+
+
+
+ template
+ void
+ BlackoilModel::computeMassFlux(const int actph ,
+ const V& transi,
+ const ADB& kr ,
+ const ADB& phasePressure,
+ const SolutionState& state)
+ {
+ const int canonicalPhaseIdx = canph_[ actph ];
+
+ const std::vector cond = phaseCondition();
+
+ const ADB tr_mult = transMult(state.pressure);
+ const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_);
+
+ rq_[ actph ].mob = tr_mult * kr / mu;
+
+ const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_);
+
+ ADB& head = rq_[ actph ].head;
+
+ // compute gravity potensial using the face average as in eclipse and MRST
+ const ADB rhoavg = ops_.caver * rho;
+
+ 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;
+
+ UpwindSelector upwind(grid_, ops_, head.value());
+
+ const ADB& b = rq_[ actph ].b;
+ const ADB& mob = rq_[ actph ].mob;
+ rq_[ actph ].mflux = upwind.select(b * mob) * head;
+ // OPM_AD_DUMP(rq_[ actph ].mob);
+ // OPM_AD_DUMP(rq_[ actph ].mflux);
+ }
+
+
+
+
+
+ template
+ void
+ BlackoilModel::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);
+
+ // Find the current sign for the threshold modification
+ const V sign_dp = sign(dp.value());
+ const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_;
+
+ // Modify potential and nullify where appropriate.
+ dp = keep_high_potential * (dp - threshold_modification);
+ }
+
+
+
+
+ template
+ std::vector
+ BlackoilModel::computeResidualNorms() const
+ {
+ std::vector residualNorms;
+
+ std::vector::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
+ const std::vector::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
+
+ for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
+ const double massBalanceResid = detail::infinityNorm( (*massBalanceIt),
+ linsolver_.parallelInformation() );
+ if (!std::isfinite(massBalanceResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residualNorms.push_back(massBalanceResid);
+ }
+
+ // the following residuals are not used in the oscillation detection now
+ const double wellFluxResid = detail::infinityNormWell( residual_.well_flux_eq,
+ linsolver_.parallelInformation() );
+ if (!std::isfinite(wellFluxResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residualNorms.push_back(wellFluxResid);
+
+ const double wellResid = detail::infinityNormWell( residual_.well_eq,
+ linsolver_.parallelInformation() );
+ if (!std::isfinite(wellResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residualNorms.push_back(wellResid);
+
+ return residualNorms;
+ }
+
+
+
+
+ template
+ double
+ BlackoilModel::convergenceReduction(const Eigen::Array& B,
+ const Eigen::Array& tempV,
+ const Eigen::Array& R,
+ std::array& R_sum,
+ std::array& maxCoeff,
+ std::array& B_avg,
+ std::vector& maxNormWell,
+ int nc,
+ int nw) const
+ {
+ // Do the global reductions
+#if HAVE_MPI
+ if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+ {
+ const ParallelISTLInformation& info =
+ boost::any_cast(linsolver_.parallelInformation());
+
+ // Compute the global number of cells and porevolume
+ std::vector v(nc, 1);
+ auto nc_and_pv = std::tuple(0, 0.0);
+ auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(),
+ Opm::Reduction::makeGlobalSumFunctor());
+ auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
+ info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
+
+ for ( int idx=0; idx(0.0 ,0.0 ,0.0);
+ auto containers = std::make_tuple(B.col(idx),
+ tempV.col(idx),
+ R.col(idx));
+ auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(),
+ Opm::Reduction::makeGlobalMaxFunctor(),
+ Opm::Reduction::makeGlobalSumFunctor());
+ info.computeReduction(containers, operators, values);
+ B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
+ maxCoeff[idx] = std::get<1>(values);
+ R_sum[idx] = std::get<2>(values);
+ maxNormWell[idx] = 0.0;
+ for ( int w=0; w(nc_and_pv);
+ }
+ else
+#endif
+ {
+ for ( int idx=0; idx
+ bool
+ BlackoilModel::getConvergence(const double dt, const int iteration)
+ {
+ const double tol_mb = param_.tolerance_mb_;
+ const double tol_cnv = param_.tolerance_cnv_;
+ const double tol_wells = param_.tolerance_wells_;
+
+ const int nc = Opm::AutoDiffGrid::numCells(grid_);
+ const int nw = wellsActive() ? wells().number_of_wells : 0;
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+ const V pv = geo_.poreVolume();
+
+ const std::vector cond = phaseCondition();
+
+ std::array CNV = {{0., 0., 0.}};
+ std::array R_sum = {{0., 0., 0.}};
+ std::array B_avg = {{0., 0., 0.}};
+ std::array maxCoeff = {{0., 0., 0.}};
+ std::array mass_balance_residual = {{0., 0., 0.}};
+ std::array well_flux_residual = {{0., 0., 0.}};
+ std::size_t cols = MaxNumPhases; // needed to pass the correct type to Eigen
+ Eigen::Array