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

Split functionality between model and solver.

Browse Source 
The step() method and everything to do with relaxation and oscillation
detection is now in the FullyImplicitSolver class.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-05-19 10:19:30 +02:00
parent 24ab95122d
commit 2e7e6c6344
5 changed files with 245 additions and 2822 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

131
opm/autodiff/BlackoilModel.hpp
View File

@ -57,6 +57,10 @@ namespace Opm {
class BlackoilModel
{
public:
// --------- Types and enums ---------
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
@ -64,28 +68,24 @@ namespace Opm {
enum RelaxType { DAMPEN, SOR };
// class holding the solver parameters
struct SolverParameter
struct ModelParameter
{
double dp_max_rel_;
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 );
SolverParameter();
ModelParameter( const parameter::ParameterGroup& param );
ModelParameter();
void reset();
};
// --------- Public methods ---------
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
@ -96,7 +96,7 @@ namespace Opm {
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
BlackoilModel(const SolverParameter& param,
BlackoilModel(const ModelParameter& param,
const Grid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
@ -105,7 +105,7 @@ namespace Opm {
const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output);
const bool terminal_output);
/// \brief Set threshold pressures that prevent or reduce flow.
/// This prevents flow across faces if the potential
@ -117,26 +117,57 @@ namespace Opm {
/// of the grid passed in the constructor.
void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
/// Take a single forward step, modifiying
/// state.pressure()
/// state.faceflux()
/// state.saturation()
/// state.gasoilratio()
/// wstate.bhp()
/// \param[in] dt time step size
/// \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);
/// Called once before each time step.
/// \param[in] dt time step size
/// \param[in] reservoir_state reservoir state variables
/// \param[in] well_state well state variables
void prepareStep(const double dt,
const ReservoirState& reservoir_state,
const WellState& well_state);
/// Assemble the residual and Jacobian of the nonlinear system.
/// \param[in] dt time step size
/// \param[in] reservoir_state reservoir state variables
/// \param[in] well_state well state variables
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;
void updateState(const V& dx,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
/// Return true if output to cout is wanted.
bool terminalOutput() const;
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
bool getConvergence(const double dt, const int iteration);
/// The number of active phases in the model.
int numPhases() const;
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
// --------- Types and enums ---------
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
@ -178,7 +209,8 @@ namespace Opm {
enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
// Member data
// --------- Data members ---------
const Grid& grid_;
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
@ -195,7 +227,7 @@ namespace Opm {
const bool has_disgas_;
const bool has_vapoil_;
SolverParameter param_;
ModelParameter param_;
bool use_threshold_pressure_;
V threshold_pressures_by_interior_face_;
@ -209,8 +241,9 @@ namespace Opm {
bool terminal_output_;
std::vector<int> primalVariable_;
V pvdt_;
// Private methods.
// --------- Private methods ---------
// return true if wells are available
bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
@ -247,18 +280,6 @@ namespace Opm {
void updateWellControls(WellStateFullyImplicitBlackoil& xw) const;
void
assemble(const V& dtpv,
const BlackoilState& x,
const bool initial_assembly,
WellStateFullyImplicitBlackoil& xw);
V solveJacobianSystem() const;
void updateState(const V& dx,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
std::vector<ADB>
computePressures(const SolutionState& state) const;
@ -286,12 +307,6 @@ namespace Opm {
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
fluidViscosity(const int phase,
const ADB& p ,
@ -365,10 +380,6 @@ namespace Opm {
void
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.
/// \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
@ -397,21 +408,9 @@ namespace Opm {
std::array<double,MaxNumPhases>& B_avg,
int nc) 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 dsMax() const { return param_.ds_max_; }
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_; }
};

281
opm/autodiff/BlackoilModel_impl.hpp
View File

@ -136,19 +136,13 @@ namespace detail {
} // namespace detail
template <class Grid>
void BlackoilModel<Grid>::SolverParameter::
void BlackoilModel<Grid>::ModelParameter::
reset()
{
// default values for the solver parameters
dp_max_rel_ = 1.0e9;
ds_max_ = 0.2;
dr_max_rel_ = 1.0e9;
relax_type_ = DAMPEN;
relax_max_ = 0.5;
relax_increment_ = 0.1;
relax_rel_tol_ = 0.2;
max_iter_ = 15; // not more then 15 its by default
min_iter_ = 1; // Default to always do at least one nonlinear iteration.
max_residual_allowed_ = 1e7;
tolerance_mb_ = 1.0e-5;
tolerance_cnv_ = 1.0e-2;
@ -156,16 +150,16 @@ namespace detail {
}
template <class Grid>
BlackoilModel<Grid>::SolverParameter::
SolverParameter()
BlackoilModel<Grid>::ModelParameter::
ModelParameter()
{
// set default values
reset();
}
template <class Grid>
BlackoilModel<Grid>::SolverParameter::
SolverParameter( const parameter::ParameterGroup& param )
BlackoilModel<Grid>::ModelParameter::
ModelParameter( const parameter::ParameterGroup& param )
{
// set default values
reset();
@ -174,38 +168,25 @@ namespace detail {
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_);
relax_max_ = param.getDefault("relax_max", relax_max_);
max_iter_ = param.getDefault("max_iter", max_iter_);
min_iter_ = param.getDefault("min_iter", min_iter_);
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_ );
std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
if (relaxation_type == "dampen") {
relax_type_ = DAMPEN;
} else if (relaxation_type == "sor") {
relax_type_ = SOR;
} else {
OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
}
}
template <class Grid>
BlackoilModel<Grid>::
BlackoilModel(const SolverParameter& 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)
BlackoilModel(const ModelParameter& 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)
@ -243,6 +224,66 @@ namespace detail {
template <class Grid>
void
BlackoilModel<Grid>::
prepareStep(const double dt,
const ReservoirState& reservoir_state,
const WellState& /* well_state */)
{
pvdt_ = geo_.poreVolume() / dt;
if (active_[Gas]) {
updatePrimalVariableFromState(reservoir_state);
}
}
template <class Grid>
int
BlackoilModel<Grid>::
sizeNonLinear() const
{
return residual_.sizeNonLinear();
}
template <class Grid>
int
BlackoilModel<Grid>::
linearIterationsLastSolve() const
{
return linsolver_.iterations();
}
template <class Grid>
bool
BlackoilModel<Grid>::
terminalOutput() const
{
return terminal_output_;
}
template <class Grid>
int
BlackoilModel<Grid>::
numPhases() const
{
return fluid_.numPhases();
}
template <class Grid>
void
BlackoilModel<Grid>::
@ -264,83 +305,6 @@ namespace detail {
template <class Grid>
int
BlackoilModel<Grid>::
step(const double dt,
BlackoilState& x ,
WellStateFullyImplicitBlackoil& xw)
{
const V pvdt = geo_.poreVolume() / dt;
if (active_[Gas]) { updatePrimalVariableFromState(x); }
// For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations
std::vector<std::vector<double>> residual_norms_history;
assemble(pvdt, x, true, xw);
bool converged = false;
double omega = 1.;
residual_norms_history.push_back(computeResidualNorms());
int it = 0;
converged = getConvergence(dt,it);
const int sizeNonLinear = residual_.sizeNonLinear();
V dxOld = V::Zero(sizeNonLinear);
bool isOscillate = false;
bool isStagnate = false;
const enum RelaxType relaxtype = relaxType();
int linearIterations = 0;
while ( (!converged && (it < maxIter())) || (minIter() > it)) {
V dx = solveJacobianSystem();
// store number of linear iterations used
linearIterations += linsolver_.iterations();
detectNewtonOscillations(residual_norms_history, it, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) {
omega -= relaxIncrement();
omega = std::max(omega, relaxMax());
if (terminal_output_)
{
std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
}
}
stablizeNewton(dx, dxOld, omega, relaxtype);
updateState(dx, x, xw);
assemble(pvdt, x, false, xw);
residual_norms_history.push_back(computeResidualNorms());
// increase iteration counter
++it;
converged = getConvergence(dt,it);
}
if (!converged) {
std::cerr << "WARNING: Failed to compute converged solution in " << it << " iterations." << std::endl;
return -1; // -1 indicates that the solver has to be restarted
}
return linearIterations;
}
template <class Grid>
BlackoilModel<Grid>::ReservoirResidualQuant::ReservoirResidualQuant()
: accum(2, ADB::null())
@ -752,19 +716,18 @@ namespace detail {
template <class Grid>
void
BlackoilModel<Grid>::
assemble(const V& pvdt,
const BlackoilState& x ,
const bool initial_assembly,
WellStateFullyImplicitBlackoil& xw )
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(xw);
updateWellControls(well_state);
// Create the primary variables.
SolutionState state = variableState(x, xw);
SolutionState state = variableState(reservoir_state, well_state);
if (initial_assembly) {
// Create the (constant, derivativeless) initial state.
@ -773,7 +736,7 @@ namespace detail {
// Compute initial accumulation contributions
// and well connection pressures.
computeAccum(state0, 0);
computeWellConnectionPressures(state0, xw);
computeWellConnectionPressures(state0, well_state);
}
// DISKVAL(state.pressure);
@ -801,18 +764,10 @@ namespace detail {
const std::vector<ADB> kr = computeRelPerm(state);
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
// std::cout << "===== kr[" << phase << "] = \n" << std::endl;
// std::cout << kr[phase];
// std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
// std::cout << rq_[phase].mflux;
residual_.material_balance_eq[ phaseIdx ] =
pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux;
// DUMP(ops_.div*rq_[phase].mflux);
// DUMP(residual_.material_balance_eq[phase]);
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@ -843,8 +798,8 @@ namespace detail {
// Add contribution from wells and set up the well equations.
V aliveWells;
addWellEq(state, xw, aliveWells);
addWellControlEq(state, xw, aliveWells);
addWellEq(state, well_state, aliveWells);
addWellControlEq(state, well_state, aliveWells);
}
@ -1756,74 +1711,6 @@ namespace detail {
return residualNorms;
}
template <class Grid>
void
BlackoilModel<Grid>::detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol,
bool& oscillate, bool& stagnate) const
{
// The detection of oscillation in two primary variable results in the report of the detection
// of oscillation for the solver.
// Only the saturations are used for oscillation detection for the black oil model.
// Stagnate is not used for any treatment here.
if ( it < 2 ) {
oscillate = false;
stagnate = false;
return;
}
stagnate = true;
int oscillatePhase = 0;
const std::vector<double>& F0 = residual_history[it];
const std::vector<double>& F1 = residual_history[it - 1];
const std::vector<double>& F2 = residual_history[it - 2];
for (int p= 0; p < fluid_.numPhases(); ++p){
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2);
// Process is 'stagnate' unless at least one phase
// exhibits significant residual change.
stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
}
oscillate = (oscillatePhase > 1);
}
template <class Grid>
void
BlackoilModel<Grid>::stablizeNewton(V& dx, V& dxOld, const double omega,
const RelaxType relax_type) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
const V tempDxOld = dxOld;
dxOld = dx;
switch (relax_type) {
case DAMPEN:
if (omega == 1.) {
return;
}
dx = dx*omega;
return;
case SOR:
if (omega == 1.) {
return;
}
dx = dx*omega + (1.-omega)*tempDxOld;
return;
default:
OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
}
return;
}
template <class Grid>
double
BlackoilModel<Grid>::convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,

401
opm/autodiff/FullyImplicitSolver.hpp
View File

@ -23,7 +23,7 @@
// #include <cassert>
// #include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
// #include <opm/autodiff/AutoDiffHelpers.hpp>
// #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
// #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
@ -36,7 +36,7 @@
namespace Opm {
// namespace parameter { class ParameterGroup; }
namespace parameter { class ParameterGroup; }
// class DerivedGeology;
// class RockCompressibility;
// class NewtonIterationBlackoilInterface;
@ -46,23 +46,52 @@ namespace Opm {
class FullyImplicitSolver
{
public:
// --------- Types and enums ---------
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
// The Newton relaxation scheme type
enum RelaxType { DAMPEN, SOR };
// Solver parameters controlling nonlinear Newton process.
struct SolverParameter
{
enum RelaxType relax_type_;
double relax_max_;
double relax_increment_;
double relax_rel_tol_;
int max_iter_; // max newton iterations
int min_iter_; // min newton iterations
SolverParameter( const parameter::ParameterGroup& param );
SolverParameter();
void reset();
};
// Forwarding types from PhysicalModel.
typedef typename PhysicalModel::ReservoirState ReservoirState;
typedef typename PhysicalModel::WellState WellState;
/// Construct solver for a given model.
explicit FullyImplicitSolver(PhysicalModel& model);
// --------- Public methods ---------
/// Take a single forward step, after which the state will be modified
/// according to PhysicalModel.
/// \param[in] dt time step size
/// \param[in] state reservoir state
/// \param[in] wstate well state
/// \return number of linear iterations used
/// Construct solver for a given model.
/// \param[in] param parameters controlling nonlinear Newton process
/// \param[in, out] model physical simulation model
explicit FullyImplicitSolver(const SolverParameter& param,
PhysicalModel& model);
/// Take a single forward step, after which the states will be modified
/// according to the physical model.
/// \param[in] dt time step size
/// \param[in] reservoir_state reservoir state variables
/// \param[in] well_state well state variables
/// \return number of linear iterations used
int
step(const double dt,
ReservoirState& state,
WellState& wstate);
ReservoirState& reservoir_state,
WellState& well_state);
/// Number of Newton iterations used in all calls to step().
unsigned int newtonIterations() const;
@ -71,351 +100,23 @@ namespace Opm {
unsigned int linearIterations() const;
private:
// --------- Data members ---------
SolverParameter param_;
PhysicalModel& model_;
unsigned int newtonIterations_;
unsigned int linearIterations_;
/*
// the Newton relaxation type
enum RelaxType { DAMPEN, SOR };
// class holding the solver parameters
struct SolverParameter
{
double dp_max_rel_;
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 );
SolverParameter();
void reset();
};
/// Construct a solver. 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
FullyImplicitSolver(const SolverParameter& 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<double>& threshold_pressures_by_face);
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
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 };
// Member data
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_;
SolverParameter 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_;
// 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;
void
assemble(const V& dtpv,
const BlackoilState& x,
const bool initial_assembly,
WellStateFullyImplicitBlackoil& xw);
V solveJacobianSystem() const;
void updateState(const V& dx,
BlackoilState& state,
WellStateFullyImplicitBlackoil& well_state);
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);
/// \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
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();
/// 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.
/// \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[in] nc The number of cells of 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,
int nc) const;
// --------- Private methods ---------
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_; }
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 dsMax() const { return param_.ds_max_; }
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_; }
*/
void stabilizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
};
} // namespace Opm

2240
opm/autodiff/FullyImplicitSolver_impl.hpp
View File

File diff suppressed because it is too large Load Diff

14
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -232,7 +232,13 @@ namespace Opm
std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
typename BlackoilModel<T>::SolverParameter solverParam( param_ );
typedef T Grid;
typedef BlackoilModel<Grid> Model;
typedef typename Model::ModelParameter ModelParam;
ModelParam modelParam( param_ );
typedef FullyImplicitSolver<Model> Solver;
typedef typename Solver::SolverParameter SolverParam;
SolverParam solverParam( param_ );
// adaptive time stepping
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
@ -292,13 +298,11 @@ namespace Opm
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
typedef T Grid;
typedef BlackoilModel<Grid> Model;
Model model(solverParam, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_);
Model model(modelParam, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_);
if (!threshold_pressures_by_face_.empty()) {
model.setThresholdPressures(threshold_pressures_by_face_);
}
FullyImplicitSolver<Model> solver(model);
Solver solver(solverParam, model);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge