OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #385 from atgeirr/polymorphism-for-blackoilmodel

Browse Source 
Static polymorphism for black-oil model
This commit is contained in:
Atgeirr Flø Rasmussen 2015-05-28 13:28:43 +02:00
commit a991eb55e3
8 changed files with 648 additions and 3283 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists_files.cmake
View File

@ -36,6 +36,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/SimulatorIncompTwophaseAd.cpp opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/BlackoilPropsAdFromDeck.cpp opm/autodiff/BlackoilPropsAdFromDeck.cpp
opm/autodiff/BlackoilModelParameters.cpp
opm/autodiff/WellDensitySegmented.cpp opm/autodiff/WellDensitySegmented.cpp
opm/autodiff/LinearisedBlackoilResidual.cpp opm/autodiff/LinearisedBlackoilResidual.cpp
) )
@ -95,7 +96,10 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AutoDiff.hpp opm/autodiff/AutoDiff.hpp
opm/autodiff/BackupRestore.hpp opm/autodiff/BackupRestore.hpp
opm/autodiff/BlackoilModel.hpp opm/autodiff/BlackoilModel.hpp
opm/autodiff/BlackoilModel_impl.hpp opm/autodiff/BlackoilModelBase.hpp
opm/autodiff/BlackoilModelBase_impl.hpp
opm/autodiff/BlackoilModelEnums.hpp
opm/autodiff/BlackoilModelParameters.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/BlackoilPropsAdInterface.hpp
opm/autodiff/CPRPreconditioner.hpp opm/autodiff/CPRPreconditioner.hpp
@ -105,8 +109,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/GeoProps.hpp opm/autodiff/GeoProps.hpp
opm/autodiff/GridHelpers.hpp opm/autodiff/GridHelpers.hpp
opm/autodiff/ImpesTPFAAD.hpp opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/FullyImplicitBlackoilSolver.hpp
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
opm/autodiff/NewtonIterationBlackoilCPR.hpp opm/autodiff/NewtonIterationBlackoilCPR.hpp
opm/autodiff/NewtonIterationBlackoilInterface.hpp opm/autodiff/NewtonIterationBlackoilInterface.hpp
opm/autodiff/NewtonIterationBlackoilSimple.hpp opm/autodiff/NewtonIterationBlackoilSimple.hpp

412
opm/autodiff/BlackoilModel.hpp
View File

@ -23,29 +23,13 @@
#ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED #ifndef OPM_BLACKOILMODEL_HEADER_INCLUDED
#define OPM_BLACKOILMODEL_HEADER_INCLUDED #define OPM_BLACKOILMODEL_HEADER_INCLUDED
#include <cassert> #include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp> #include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <array>
struct UnstructuredGrid;
struct Wells;
namespace Opm { namespace Opm {
namespace parameter { class ParameterGroup; }
class DerivedGeology;
class RockCompressibility;
class NewtonIterationBlackoilInterface;
class BlackoilState;
class WellStateFullyImplicitBlackoil;
/// A model implementation for three-phase black oil. /// A model implementation for three-phase black oil.
/// ///
/// The simulator is capable of handling three-phase problems /// The simulator is capable of handling three-phase problems
@ -56,34 +40,10 @@ namespace Opm {
/// It uses automatic differentiation via the class AutoDiffBlock /// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix. /// to simplify assembly of the jacobian matrix.
template<class Grid> template<class Grid>
class BlackoilModel class BlackoilModel : public BlackoilModelBase<Grid, BlackoilModel<Grid> >
{ {
public: public:
// --------- Types and enums --------- typedef BlackoilModelBase<Grid, BlackoilModel<Grid> > Base;
typedef AutoDiffBlock<double> 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 /// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to /// arguments of this functions, and they are expected to
@ -98,344 +58,34 @@ namespace Opm {
/// \param[in] has_disgas turn on dissolved gas /// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature /// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr /// \param[in] terminal_output request output to cout/cerr
BlackoilModel(const ModelParameters& param, BlackoilModel(const typename Base::ModelParameters& param,
const Grid& grid , const Grid& grid,
const BlackoilPropsAdInterface& fluid, const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo , const DerivedGeology& geo,
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
const Wells* wells, const Wells* wells,
const NewtonIterationBlackoilInterface& linsolver, const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas, const bool has_disgas,
const bool has_vapoil, const bool has_vapoil,
const bool terminal_output); const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
/// \brief Set threshold pressures that prevent or reduce flow. has_disgas, has_vapoil, terminal_output)
/// 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);
/// 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<double> 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<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
struct ReservoirResidualQuant {
ReservoirResidualQuant();
std::vector<ADB> 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<ADB> saturation;
ADB rs;
ADB rv;
ADB qs;
ADB bhp;
// Below are quantities stored in the state for optimization purposes.
std::vector<ADB> 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<bool> active_;
// Size = # active phases. Maps active -> canonical phase indices.
const std::vector<int> canph_;
const std::vector<int> 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<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> 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<int> 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<ADB>
computePressures(const SolutionState& state) const;
std::vector<ADB>
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<ADB>
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<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const;
V
fluidRsSat(const V& p,
const V& so,
const std::vector<int>& cells) const;
ADB
fluidRsSat(const ADB& p,
const ADB& so,
const std::vector<int>& cells) const;
V
fluidRvSat(const V& p,
const V& so,
const std::vector<int>& cells) const;
ADB
fluidRvSat(const ADB& p,
const ADB& so,
const std::vector<int>& cells) const;
ADB
poroMult(const ADB& p) const;
ADB
transMult(const ADB& p) const;
void
classifyCondition(const SolutionState& state,
std::vector<PhasePresence>& cond ) const;
const std::vector<PhasePresence>
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<double, Eigen::Dynamic, MaxNumPhases>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& R,
std::array<double,MaxNumPhases>& R_sum,
std::array<double,MaxNumPhases>& maxCoeff,
std::array<double,MaxNumPhases>& B_avg,
std::vector<double>& 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_; }
}; };
/// Providing types by template specialisation of ModelTraits for BlackoilModel.
template <class Grid>
struct ModelTraits< BlackoilModel<Grid> >
{
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
};
} // namespace Opm } // namespace Opm
#include "BlackoilModel_impl.hpp"
#endif // OPM_BLACKOILMODEL_HEADER_INCLUDED #endif // OPM_BLACKOILMODEL_HEADER_INCLUDED

334
opm/autodiff/FullyImplicitBlackoilSolver.hpp → opm/autodiff/BlackoilModelBase.hpp
View File

@ -1,7 +1,7 @@
/* /*
Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014, 2015 Statoil ASA.
Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
Copyright 2014, 2015 Statoil AS
Copyright 2015 NTNU Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -20,8 +20,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED #ifndef OPM_BLACKOILMODELBASE_HEADER_INCLUDED
#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED #define OPM_BLACKOILMODELBASE_HEADER_INCLUDED
#include <cassert> #include <cassert>
@ -30,10 +30,10 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp> #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp> #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <array> #include <array>
struct UnstructuredGrid;
struct Wells; struct Wells;
namespace Opm { namespace Opm {
@ -42,52 +42,72 @@ namespace Opm {
class DerivedGeology; class DerivedGeology;
class RockCompressibility; class RockCompressibility;
class NewtonIterationBlackoilInterface; class NewtonIterationBlackoilInterface;
class BlackoilState;
class WellStateFullyImplicitBlackoil;
/// A fully implicit solver for the black-oil problem. /// Struct for containing iteration variables.
struct DefaultBlackoilSolutionState
{
typedef AutoDiffBlock<double> ADB;
explicit DefaultBlackoilSolutionState(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())
{
}
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
ADB qs;
ADB bhp;
// Below are quantities stored in the state for optimization purposes.
std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
};
/// Traits to encapsulate the types used by classes using or
/// extending this model. Forward declared here, must be
/// specialised for each concrete model class.
template <class ConcreteModel>
struct ModelTraits;
/// A model implementation for three-phase black oil.
/// ///
/// The simulator is capable of handling three-phase problems /// The simulator is capable of handling three-phase problems
/// where gas can be dissolved in oil (but not vice versa). It /// where gas can be dissolved in oil and vice versa. It
/// uses an industry-standard TPFA discretization with per-phase /// uses an industry-standard TPFA discretization with per-phase
/// upwind weighting of mobilities. /// upwind weighting of mobilities.
/// ///
/// It uses automatic differentiation via the class AutoDiffBlock /// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix. /// to simplify assembly of the jacobian matrix.
template<class T> /// \tparam Grid UnstructuredGrid or CpGrid.
class FullyImplicitBlackoilSolver /// \tparam Implementation Provides concrete state types.
template<class Grid, class Implementation>
class BlackoilModelBase
{ {
public: public:
// the Newton relaxation type // --------- Types and enums ---------
enum RelaxType { DAMPEN, SOR }; typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
// class holding the solver parameters typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;
struct SolverParameter typedef typename ModelTraits<Implementation>::WellState WellState;
{ typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
double dp_max_rel_; typedef typename ModelTraits<Implementation>::SolutionState SolutionState;
double ds_max_;
double dr_max_rel_;
enum RelaxType relax_type_;
double relax_max_;
double relax_increment_;
double relax_rel_tol_;
double max_residual_allowed_;
double tolerance_mb_;
double tolerance_cnv_;
double tolerance_wells_;
int max_iter_; // max newton iterations
int min_iter_; // min newton iterations
SolverParameter( const parameter::ParameterGroup& param ); // --------- Public methods ---------
SolverParameter();
void reset(); /// Construct the model. It will retain references to the
};
/// \brief The type of the grid that we use.
typedef T Grid;
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to /// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver. /// remain in scope for the lifetime of the solver.
/// \param[in] param parameters /// \param[in] param parameters
@ -97,16 +117,19 @@ namespace Opm {
/// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure /// \param[in] wells well structure
/// \param[in] linsolver linear solver /// \param[in] linsolver linear solver
FullyImplicitBlackoilSolver(const SolverParameter& param, /// \param[in] has_disgas turn on dissolved gas
const Grid& grid , /// \param[in] has_vapoil turn on vaporized oil feature
const BlackoilPropsAdInterface& fluid, /// \param[in] terminal_output request output to cout/cerr
const DerivedGeology& geo , BlackoilModelBase(const ModelParameters& param,
const RockCompressibility* rock_comp_props, const Grid& grid ,
const Wells* wells, const BlackoilPropsAdInterface& fluid,
const NewtonIterationBlackoilInterface& linsolver, const DerivedGeology& geo ,
const bool has_disgas, const RockCompressibility* rock_comp_props,
const bool has_vapoil, const Wells* wells,
const bool terminal_output); const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output);
/// \brief Set threshold pressures that prevent or reduce flow. /// \brief Set threshold pressures that prevent or reduce flow.
/// This prevents flow across faces if the potential /// This prevents flow across faces if the potential
@ -118,29 +141,71 @@ namespace Opm {
/// of the grid passed in the constructor. /// of the grid passed in the constructor.
void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face); void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
/// Take a single forward step, modifiying /// Called once before each time step.
/// state.pressure() /// \param[in] dt time step size
/// state.faceflux() /// \param[in, out] reservoir_state reservoir state variables
/// state.saturation() /// \param[in, out] well_state well state variables
/// state.gasoilratio() void prepareStep(const double dt,
/// wstate.bhp() ReservoirState& reservoir_state,
/// \param[in] dt time step size WellState& well_state);
/// \param[in] state reservoir state
/// \param[in] wstate well state
/// \return number of linear iterations used
int
step(const double dt ,
BlackoilState& state ,
WellStateFullyImplicitBlackoil& wstate);
unsigned int newtonIterations () const { return newtonIterations_; } /// Called once after each time step.
unsigned int linearIterations () const { return linearIterations_; } /// 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 ReservoirState& reservoir_state,
WellState& 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<double> 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,
ReservoirState& reservoir_state,
WellState& 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;
protected:
// --------- Types and enums ---------
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double, typedef Eigen::Array<double,
Eigen::Dynamic, Eigen::Dynamic,
Eigen::Dynamic, Eigen::Dynamic,
@ -151,38 +216,18 @@ namespace Opm {
std::vector<ADB> accum; // Accumulations std::vector<ADB> accum; // Accumulations
ADB mflux; // Mass flux (surface conditions) ADB mflux; // Mass flux (surface conditions)
ADB b; // Reciprocal FVF ADB b; // Reciprocal FVF
ADB head; // Pressure drop across int. interfaces ADB dh; // Pressure drop across int. interfaces
ADB mob; // Phase mobility (per cell) ADB mob; // Phase mobility (per cell)
}; };
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
ADB qs;
ADB bhp;
// Below are quantities stored in the state for optimization purposes.
std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
};
struct WellOps { struct WellOps {
WellOps(const Wells* wells); WellOps(const Wells* wells);
M w2p; // well -> perf (scatter) M w2p; // well -> perf (scatter)
M p2w; // perf -> well (gather) M p2w; // perf -> well (gather)
}; };
enum { Water = BlackoilPropsAdInterface::Water, // --------- Data members ---------
Oil = BlackoilPropsAdInterface::Oil ,
Gas = BlackoilPropsAdInterface::Gas ,
MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
};
enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
// Member data
const Grid& grid_; const Grid& grid_;
const BlackoilPropsAdInterface& fluid_; const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_; const DerivedGeology& geo_;
@ -199,74 +244,94 @@ namespace Opm {
const bool has_disgas_; const bool has_disgas_;
const bool has_vapoil_; const bool has_vapoil_;
SolverParameter param_; ModelParameters param_;
bool use_threshold_pressure_; bool use_threshold_pressure_;
V threshold_pressures_by_interior_face_; V threshold_pressures_by_interior_face_;
std::vector<ReservoirResidualQuant> rq_; std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_; std::vector<PhasePresence> phaseCondition_;
V isRs_;
V isRv_;
V isSg_;
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
LinearisedBlackoilResidual residual_; LinearisedBlackoilResidual residual_;
/// \brief Whether we print something to std::cout /// \brief Whether we print something to std::cout
bool terminal_output_; bool terminal_output_;
unsigned int newtonIterations_;
unsigned int linearIterations_;
std::vector<int> primalVariable_; std::vector<int> primalVariable_;
V pvdt_;
// Private methods. // --------- Protected methods ---------
/// Access the most-derived class used for
/// static polymorphism (CRTP).
Implementation& asImpl()
{
return static_cast<Implementation&>(*this);
}
/// Access the most-derived class used for
/// static polymorphism (CRTP).
const Implementation& asImpl() const
{
return static_cast<const Implementation&>(*this);
}
// return true if wells are available // return true if wells are available
bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; } bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
// return wells object // return wells object
const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; } const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
SolutionState
constantState(const BlackoilState& x,
const WellStateFullyImplicitBlackoil& xw) const;
void void
makeConstantState(SolutionState& state) const; makeConstantState(SolutionState& state) const;
SolutionState SolutionState
variableState(const BlackoilState& x, variableState(const ReservoirState& x,
const WellStateFullyImplicitBlackoil& xw) const; const WellState& xw) const;
std::vector<V>
variableStateInitials(const ReservoirState& x,
const WellState& xw) const;
std::vector<int>
variableStateIndices() const;
SolutionState
variableStateExtractVars(const ReservoirState& x,
const std::vector<int>& indices,
std::vector<ADB>& vars) const;
void void
computeAccum(const SolutionState& state, computeAccum(const SolutionState& state,
const int aix ); const int aix );
void computeWellConnectionPressures(const SolutionState& state, void computeWellConnectionPressures(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw); const WellState& xw);
void
assembleMassBalanceEq(const SolutionState& state);
void void
addWellControlEq(const SolutionState& state, addWellControlEq(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw, const WellState& xw,
const V& aliveWells); const V& aliveWells);
void void
addWellEq(const SolutionState& state, addWellEq(const SolutionState& state,
WellStateFullyImplicitBlackoil& xw, WellState& xw,
V& aliveWells); V& aliveWells);
void updateWellControls(WellStateFullyImplicitBlackoil& xw) const;
void void
assemble(const V& dtpv, extraAddWellEq(const SolutionState& state,
const BlackoilState& x, const WellState& xw,
const bool initial_assembly, const std::vector<ADB>& cq_ps,
WellStateFullyImplicitBlackoil& xw); const std::vector<ADB>& cmix_s,
const ADB& cqt_is,
const std::vector<int>& well_cells);
V solveJacobianSystem() const; void updateWellControls(WellState& xw) const;
void updateState(const V& dx,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
std::vector<ADB>
computePressures(const SolutionState& state) const;
std::vector<ADB> std::vector<ADB>
computePressures(const ADB& po, computePressures(const ADB& po,
@ -292,12 +357,6 @@ namespace Opm {
void applyThresholdPressures(ADB& dp); void applyThresholdPressures(ADB& dp);
/// \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<double> computeResidualNorms() const;
ADB ADB
fluidViscosity(const int phase, fluidViscosity(const int phase,
const ADB& p , const ADB& p ,
@ -351,30 +410,23 @@ namespace Opm {
ADB ADB
transMult(const ADB& p) const; transMult(const ADB& p) const;
void
classifyCondition(const SolutionState& state,
std::vector<PhasePresence>& cond ) const;
const std::vector<PhasePresence> const std::vector<PhasePresence>
phaseCondition() const {return phaseCondition_;} phaseCondition() const {return phaseCondition_;}
void void
classifyCondition(const BlackoilState& state); classifyCondition(const ReservoirState& state);
/// update the primal variable for Sg, Rv or Rs. The Gas phase must /// update the primal variable for Sg, Rv or Rs. The Gas phase must
/// be active to call this method. /// be active to call this method.
void void
updatePrimalVariableFromState(const BlackoilState& state); updatePrimalVariableFromState(const ReservoirState& state);
/// Update the phaseCondition_ member based on the primalVariable_ member. /// Update the phaseCondition_ member based on the primalVariable_ member.
/// Also updates isRs_, isRv_ and isSg_;
void void
updatePhaseCondFromPrimalVariable(); updatePhaseCondFromPrimalVariable();
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
bool getConvergence(const double dt, const int iteration);
/// \brief Compute the reduction within the convergence check. /// \brief Compute the reduction within the convergence check.
/// \param[in] B A matrix with MaxNumPhases columns and the same number rows /// \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 /// as the number of cells of the grid. B.col(i) contains the values
@ -407,26 +459,14 @@ namespace Opm {
int nc, int nc,
int nw) const; int nw) const;
void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol,
bool& oscillate, bool& stagnate) const;
void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
double dpMaxRel() const { return param_.dp_max_rel_; } double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; } double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; } double drMaxRel() const { return param_.dr_max_rel_; }
enum RelaxType relaxType() const { return param_.relax_type_; }
double relaxMax() const { return param_.relax_max_; };
double relaxIncrement() const { return param_.relax_increment_; };
double relaxRelTol() const { return param_.relax_rel_tol_; };
double maxIter() const { return param_.max_iter_; }
double minIter() const { return param_.min_iter_; }
double maxResidualAllowed() const { return param_.max_residual_allowed_; } double maxResidualAllowed() const { return param_.max_residual_allowed_; }
}; };
} // namespace Opm } // namespace Opm
#include "FullyImplicitBlackoilSolver_impl.hpp" #include "BlackoilModelBase_impl.hpp"
#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED #endif // OPM_BLACKOILMODELBASE_HEADER_INCLUDED

601
opm/autodiff/BlackoilModel_impl.hpp → opm/autodiff/BlackoilModelBase_impl.hpp
View File

File diff suppressed because it is too large Load Diff

53
opm/autodiff/BlackoilModelEnums.hpp Normal file
View File

@ -0,0 +1,53 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILMODELENUMS_HEADER_INCLUDED
#define OPM_BLACKOILMODELENUMS_HEADER_INCLUDED
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
namespace Opm
{
enum Phases {
Water = BlackoilPropsAdInterface::Water,
Oil = BlackoilPropsAdInterface::Oil ,
Gas = BlackoilPropsAdInterface::Gas ,
MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
};
enum PrimalVariables {
Sg = 0,
RS = 1,
RV = 2
};
enum CanonicalVariablePositions {
Pressure = 0,
Sw = 1,
Xvar = 2,
Qs = 3,
Bhp = 4,
Next // For extension.
};
} // namespace Opm
#endif // OPM_BLACKOILMODELENUMS_HEADER_INCLUDED

67
opm/autodiff/BlackoilModelParameters.cpp Normal file
View File

@ -0,0 +1,67 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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 <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm
{
BlackoilModelParameters::BlackoilModelParameters()
{
// set default values
reset();
}
BlackoilModelParameters::BlackoilModelParameters( 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_ );
}
void BlackoilModelParameters::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;
}
} // namespace Opm

58
opm/autodiff/BlackoilModelParameters.hpp Normal file
View File

@ -0,0 +1,58 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED
#define OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED
namespace Opm
{
namespace parameter { class ParameterGroup; }
/// Solver parameters for the BlackoilModel.
struct BlackoilModelParameters
{
/// Max relative change in pressure in single iteration.
double dp_max_rel_;
/// Max absolute change in saturation in single iteration.
double ds_max_;
/// Max relative change in gas-oil or oil-gas ratio in single iteration.
double dr_max_rel_;
/// Absolute max limit for residuals.
double max_residual_allowed_;
/// Relative mass balance tolerance (total mass balance error).
double tolerance_mb_;
/// Local convergence tolerance (max of local saturation errors).
double tolerance_cnv_;
/// Well convergence tolerance.
double tolerance_wells_;
/// Construct from user parameters or defaults.
explicit BlackoilModelParameters( const parameter::ParameterGroup& param );
/// Construct with default parameters.
BlackoilModelParameters();
/// Set default parameters.
void reset();
};
} // namespace Opm
#endif // OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED

2398
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

File diff suppressed because it is too large Load Diff