From 620ef2a3dde9564ed684746d42e9351aa2471a15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 24 May 2015 09:59:40 +0200 Subject: [PATCH] Create BlackoilModelBase class. The class is identical to BlackoilModel class at this stage, but since it was renamed from FullyImplicitBlackoilSolver it keeps the commit history better. --- CMakeLists_files.cmake | 4 +- ...ackoilSolver.hpp => BlackoilModelBase.hpp} | 211 +++---- ...er_impl.hpp => BlackoilModelBase_impl.hpp} | 563 ++++++++---------- 3 files changed, 349 insertions(+), 429 deletions(-) rename opm/autodiff/{FullyImplicitBlackoilSolver.hpp => BlackoilModelBase.hpp} (75%) rename opm/autodiff/{FullyImplicitBlackoilSolver_impl.hpp => BlackoilModelBase_impl.hpp} (84%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5be7b04b9..44f6e33fc 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -96,6 +96,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BackupRestore.hpp opm/autodiff/BlackoilModel.hpp opm/autodiff/BlackoilModel_impl.hpp + opm/autodiff/BlackoilModelBase.hpp + opm/autodiff/BlackoilModelBase_impl.hpp opm/autodiff/BlackoilPropsAdFromDeck.hpp opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/CPRPreconditioner.hpp @@ -105,8 +107,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/GeoProps.hpp opm/autodiff/GridHelpers.hpp opm/autodiff/ImpesTPFAAD.hpp - opm/autodiff/FullyImplicitBlackoilSolver.hpp - opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp opm/autodiff/NewtonIterationBlackoilCPR.hpp opm/autodiff/NewtonIterationBlackoilInterface.hpp opm/autodiff/NewtonIterationBlackoilSimple.hpp diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/BlackoilModelBase.hpp similarity index 75% rename from opm/autodiff/FullyImplicitBlackoilSolver.hpp rename to opm/autodiff/BlackoilModelBase.hpp index f67e38b30..46ea99683 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -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 Statoil AS Copyright 2015 NTNU This file is part of the Open Porous Media project (OPM). @@ -20,8 +20,8 @@ along with OPM. If not, see . */ -#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED -#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED +#ifndef OPM_BLACKOILMODELBASE_HEADER_INCLUDED +#define OPM_BLACKOILMODELBASE_HEADER_INCLUDED #include @@ -46,48 +46,46 @@ namespace Opm { class WellStateFullyImplicitBlackoil; - /// A fully implicit solver for the black-oil problem. + /// A model implementation for three-phase black oil. /// /// 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 /// upwind weighting of mobilities. /// /// It uses automatic differentiation via the class AutoDiffBlock /// to simplify assembly of the jacobian matrix. - template - class FullyImplicitBlackoilSolver + template + class BlackoilModelBase { public: - // the Newton relaxation type - enum RelaxType { DAMPEN, SOR }; + // --------- Types and enums --------- + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + typedef BlackoilState ReservoirState; + typedef WellStateFullyImplicitBlackoil WellState; - // class holding the solver parameters - struct SolverParameter + /// Model-specific solver parameters. + struct ModelParameters { - 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 + double dp_max_rel_; + double ds_max_; + double dr_max_rel_; + double max_residual_allowed_; + double tolerance_mb_; + double tolerance_cnv_; + double tolerance_wells_; - SolverParameter( const parameter::ParameterGroup& param ); - SolverParameter(); + explicit ModelParameters( const parameter::ParameterGroup& param ); + ModelParameters(); void reset(); }; - /// \brief The type of the grid that we use. - typedef T Grid; - /// Construct a solver. It will retain references to the + // --------- 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 @@ -97,16 +95,19 @@ namespace Opm { /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] wells well structure /// \param[in] linsolver linear solver - FullyImplicitBlackoilSolver(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); + /// \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 + BlackoilModelBase(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 @@ -118,29 +119,71 @@ namespace Opm { /// of the grid passed in the constructor. void setThresholdPressures(const std::vector& threshold_pressures_by_face); - /// Take a single forward step, modifiying - /// state.pressure() - /// state.faceflux() - /// 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, 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); - unsigned int newtonIterations () const { return newtonIterations_; } - unsigned int linearIterations () const { return linearIterations_; } + /// 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 AutoDiffBlock ADB; - typedef ADB::V V; - typedef ADB::M M; + + // --------- Types and enums --------- + typedef Eigen::Array primalVariable_; + V pvdt_; - // Private methods. + // --------- Private methods --------- // return true if wells are available bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; } @@ -253,18 +296,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 computePressures(const SolutionState& state) const; @@ -292,12 +323,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 computeResidualNorms() const; - ADB fluidViscosity(const int phase, const ADB& p , @@ -371,10 +396,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 @@ -407,26 +428,14 @@ namespace Opm { int nc, int nw) const; - void detectNewtonOscillations(const std::vector>& 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_; } }; } // namespace Opm -#include "FullyImplicitBlackoilSolver_impl.hpp" +#include "BlackoilModelBase_impl.hpp" -#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED +#endif // OPM_BLACKOILMODELBASE_HEADER_INCLUDED diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp similarity index 84% rename from opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp rename to opm/autodiff/BlackoilModelBase_impl.hpp index 5362d142a..19ad7b421 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1,7 +1,7 @@ /* - Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2014, 2015 Statoil AS + Copyright 2014, 2015 Statoil ASA. Copyright 2015 NTNU Copyright 2015 IRIS AS @@ -21,7 +21,10 @@ along with OPM. If not, see . */ -#include +#ifndef OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED +#define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED + +#include #include #include @@ -50,21 +53,21 @@ //#include // A debugging utility. -#define DUMP(foo) \ +#define OPM_AD_DUMP(foo) \ do { \ std::cout << "==========================================\n" \ << #foo ":\n" \ << collapseJacs(foo) << std::endl; \ } while (0) -#define DUMPVAL(foo) \ +#define OPM_AD_DUMPVAL(foo) \ do { \ std::cout << "==========================================\n" \ << #foo ":\n" \ << foo.value() << std::endl; \ } while (0) -#define DISKVAL(foo) \ +#define OPM_AD_DISKVAL(foo) \ do { \ std::ofstream os(#foo); \ os.precision(16); \ @@ -133,37 +136,31 @@ namespace detail { } // namespace detail - template - void FullyImplicitBlackoilSolver::SolverParameter:: + template + void BlackoilModelBase::ModelParameters:: 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; tolerance_wells_ = 5.0e-1; } - template - FullyImplicitBlackoilSolver::SolverParameter:: - SolverParameter() + template + BlackoilModelBase::ModelParameters:: + ModelParameters() { // set default values reset(); } - template - FullyImplicitBlackoilSolver::SolverParameter:: - SolverParameter( const parameter::ParameterGroup& param ) + template + BlackoilModelBase::ModelParameters:: + ModelParameters( const parameter::ParameterGroup& param ) { // set default values reset(); @@ -172,38 +169,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 - FullyImplicitBlackoilSolver:: - FullyImplicitBlackoilSolver(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) + template + BlackoilModelBase:: + BlackoilModelBase(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) @@ -225,8 +209,6 @@ namespace detail { ADB::null(), ADB::null() } ) , terminal_output_ (terminal_output) - , newtonIterations_( 0 ) - , linearIterations_( 0 ) { #if HAVE_MPI if ( terminal_output_ ) { @@ -243,9 +225,82 @@ namespace detail { - template + template void - FullyImplicitBlackoilSolver:: + BlackoilModelBase:: + prepareStep(const double dt, + ReservoirState& reservoir_state, + WellState& /* well_state */) + { + pvdt_ = geo_.poreVolume() / dt; + if (active_[Gas]) { + updatePrimalVariableFromState(reservoir_state); + } + } + + + + + template + void + BlackoilModelBase:: + afterStep(const double /* dt */, + ReservoirState& /* reservoir_state */, + WellState& /* well_state */) + { + // Does nothing in this model. + } + + + + + template + int + BlackoilModelBase:: + sizeNonLinear() const + { + return residual_.sizeNonLinear(); + } + + + + + template + int + BlackoilModelBase:: + linearIterationsLastSolve() const + { + return linsolver_.iterations(); + } + + + + + template + bool + BlackoilModelBase:: + terminalOutputEnabled() const + { + return terminal_output_; + } + + + + + template + int + BlackoilModelBase:: + numPhases() const + { + return fluid_.numPhases(); + } + + + + + template + void + BlackoilModelBase:: setThresholdPressures(const std::vector& threshold_pressures_by_face) { const int num_faces = AutoDiffGrid::numFaces(grid_); @@ -264,88 +319,8 @@ namespace detail { - template - int - FullyImplicitBlackoilSolver:: - 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> 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 - } - - linearIterations_ += linearIterations; - newtonIterations_ += it; - - return linearIterations; - } - - - - - - template - FullyImplicitBlackoilSolver::ReservoirResidualQuant::ReservoirResidualQuant() + template + BlackoilModelBase::ReservoirResidualQuant::ReservoirResidualQuant() : accum(2, ADB::null()) , mflux( ADB::null()) , b ( ADB::null()) @@ -358,8 +333,8 @@ namespace detail { - template - FullyImplicitBlackoilSolver::SolutionState::SolutionState(const int np) + template + BlackoilModelBase::SolutionState::SolutionState(const int np) : pressure ( ADB::null()) , temperature( ADB::null()) , saturation(np, ADB::null()) @@ -375,8 +350,8 @@ namespace detail { - template - FullyImplicitBlackoilSolver:: + template + BlackoilModelBase:: WellOps::WellOps(const Wells* wells) : w2p(), p2w() @@ -411,9 +386,9 @@ namespace detail { - template - typename FullyImplicitBlackoilSolver::SolutionState - FullyImplicitBlackoilSolver::constantState(const BlackoilState& x, + template + typename BlackoilModelBase::SolutionState + BlackoilModelBase::constantState(const BlackoilState& x, const WellStateFullyImplicitBlackoil& xw) const { auto state = variableState(x, xw); @@ -425,9 +400,9 @@ namespace detail { - template + template void - FullyImplicitBlackoilSolver::makeConstantState(SolutionState& state) const + BlackoilModelBase::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 @@ -454,9 +429,9 @@ namespace detail { - template - typename FullyImplicitBlackoilSolver::SolutionState - FullyImplicitBlackoilSolver::variableState(const BlackoilState& x, + template + typename BlackoilModelBase::SolutionState + BlackoilModelBase::variableState(const BlackoilState& x, const WellStateFullyImplicitBlackoil& xw) const { using namespace Opm::AutoDiffGrid; @@ -612,9 +587,9 @@ namespace detail { - template + template void - FullyImplicitBlackoilSolver::computeAccum(const SolutionState& state, + BlackoilModelBase::computeAccum(const SolutionState& state, const int aix ) { const Opm::PhaseUsage& pu = fluid_.phaseUsage(); @@ -635,8 +610,8 @@ namespace detail { 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]; - // DUMP(rq_[pos].b); - // DUMP(rq_[pos].accum[aix]); + // OPM_AD_DUMP(rq_[pos].b); + // OPM_AD_DUMP(rq_[pos].accum[aix]); } } @@ -651,7 +626,7 @@ namespace detail { rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; rq_[po].accum[aix] += state.rv * accum_gas_copy; - //DUMP(rq_[pg].accum[aix]); + // OPM_AD_DUMP(rq_[pg].accum[aix]); } } @@ -659,8 +634,8 @@ namespace detail { - template - void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state, + template + void BlackoilModelBase::computeWellConnectionPressures(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw) { if( ! wellsActive() ) return ; @@ -752,22 +727,21 @@ namespace detail { - template + template void - FullyImplicitBlackoilSolver:: - assemble(const V& pvdt, - const BlackoilState& x , - const bool initial_assembly, - WellStateFullyImplicitBlackoil& xw ) + BlackoilModelBase:: + 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. @@ -776,17 +750,17 @@ namespace detail { // Compute initial accumulation contributions // and well connection pressures. computeAccum(state0, 0); - computeWellConnectionPressures(state0, xw); + computeWellConnectionPressures(state0, well_state); } - // DISKVAL(state.pressure); - // DISKVAL(state.saturation[0]); - // DISKVAL(state.saturation[1]); - // DISKVAL(state.saturation[2]); - // DISKVAL(state.rs); - // DISKVAL(state.rv); - // DISKVAL(state.qs); - // DISKVAL(state.bhp); + // 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 -------- @@ -804,18 +778,10 @@ namespace detail { 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); - // 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 -------- @@ -838,7 +804,7 @@ namespace detail { residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux); residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux); - // DUMP(residual_.material_balance_eq[ Gas ]); + // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]); } @@ -846,16 +812,16 @@ 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); } - template - void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state, + template + void BlackoilModelBase::addWellEq(const SolutionState& state, WellStateFullyImplicitBlackoil& xw, V& aliveWells) { @@ -1096,8 +1062,8 @@ namespace detail { - template - void FullyImplicitBlackoilSolver::updateWellControls(WellStateFullyImplicitBlackoil& xw) const + template + void BlackoilModelBase::updateWellControls(WellStateFullyImplicitBlackoil& xw) const { if( ! wellsActive() ) return ; @@ -1172,8 +1138,8 @@ namespace detail { - template - void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state, + template + void BlackoilModelBase::addWellControlEq(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw, const V& aliveWells) { @@ -1236,15 +1202,15 @@ namespace detail { } Selector alive_selector(aliveWells, Selector::NotEqualZero); residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs); - // DUMP(residual_.well_eq); + // OPM_AD_DUMP(residual_.well_eq); } - template - V FullyImplicitBlackoilSolver::solveJacobianSystem() const + template + V BlackoilModelBase::solveJacobianSystem() const { return linsolver_.computeNewtonIncrement(residual_); } @@ -1260,7 +1226,7 @@ namespace detail { /// \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() ) + double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() ) { #if HAVE_MPI if ( pinfo.type() == typeid(ParallelISTLInformation) ) @@ -1309,10 +1275,10 @@ namespace detail { - template - void FullyImplicitBlackoilSolver::updateState(const V& dx, - BlackoilState& state, - WellStateFullyImplicitBlackoil& well_state) + template + void BlackoilModelBase::updateState(const V& dx, + BlackoilState& reservoir_state, + WellStateFullyImplicitBlackoil& well_state) { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); @@ -1361,16 +1327,16 @@ namespace detail { // Pressure update. const double dpmaxrel = dpMaxRel(); - const V p_old = Eigen::Map(&state.pressure()[0], nc, 1); + 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, state.pressure().begin()); + 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(& state.saturation()[0], nc, np); + const DataBlock s_old = Eigen::Map(& reservoir_state.saturation()[0], nc, np); const double dsmax = dsMax(); V so; @@ -1448,19 +1414,19 @@ namespace detail { so = so / sumSat; sg = sg / sumSat; - // Update the state + // Update the reservoir_state for (int c = 0; c < nc; ++c) { - state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c]; + reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c]; } for (int c = 0; c < nc; ++c) { - state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[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) { - state.saturation()[c*np + pos] = so[c]; + reservoir_state.saturation()[c*np + pos] = so[c]; } } @@ -1468,14 +1434,14 @@ namespace detail { const double drmaxrel = drMaxRel(); V rs; if (has_disgas_) { - const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); + 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(&state.rv()[0], nc); + 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; @@ -1496,7 +1462,7 @@ namespace detail { auto hasGas = (sg > 0 && isRs == 0); // Set oil saturated if previous rs is sufficiently large - const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); + 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) { @@ -1522,7 +1488,7 @@ namespace detail { auto hasOil = (so > 0 && isRv == 0); // Set oil saturated if previous rv is sufficiently large - const V rv_old = Eigen::Map(&state.rv()[0], nc); + 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) { @@ -1535,13 +1501,13 @@ namespace detail { } - // Update the state + // Update the reservoir_state if (has_disgas_) { - std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); + std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin()); } if (has_vapoil_) { - std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); + std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin()); } if( wellsActive() ) @@ -1571,9 +1537,9 @@ namespace detail { - template + template std::vector - FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state) const + BlackoilModelBase::computeRelPerm(const SolutionState& state) const { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -1600,9 +1566,9 @@ namespace detail { - template + template std::vector - FullyImplicitBlackoilSolver::computePressures(const SolutionState& state) const + BlackoilModelBase::computePressures(const SolutionState& state) const { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -1630,9 +1596,9 @@ namespace detail { - template + template std::vector - FullyImplicitBlackoilSolver:: + BlackoilModelBase:: computePressures(const ADB& po, const ADB& sw, const ADB& so, @@ -1667,9 +1633,9 @@ namespace detail { - template + template V - FullyImplicitBlackoilSolver::computeGasPressure(const V& po, + BlackoilModelBase::computeGasPressure(const V& po, const V& sw, const V& so, const V& sg) const @@ -1686,9 +1652,9 @@ namespace detail { - template + template void - FullyImplicitBlackoilSolver::computeMassFlux(const int actph , + BlackoilModelBase::computeMassFlux(const int actph , const V& transi, const ADB& kr , const ADB& phasePressure, @@ -1724,17 +1690,17 @@ namespace detail { const ADB& b = rq_[ actph ].b; const ADB& mob = rq_[ actph ].mob; rq_[ actph ].mflux = upwind.select(b * mob) * head; - // DUMP(rq_[ actph ].mob); - // DUMP(rq_[ actph ].mflux); + // OPM_AD_DUMP(rq_[ actph ].mob); + // OPM_AD_DUMP(rq_[ actph ].mflux); } - template + template void - FullyImplicitBlackoilSolver::applyThresholdPressures(ADB& dp) + BlackoilModelBase::applyThresholdPressures(ADB& dp) { // We support reversible threshold pressures only. // Method: if the potential difference is lower (in absolute @@ -1760,9 +1726,12 @@ namespace detail { dp = keep_high_potential * (dp - threshold_modification); } - template + + + + template std::vector - FullyImplicitBlackoilSolver::computeResidualNorms() const + BlackoilModelBase::computeResidualNorms() const { std::vector residualNorms; @@ -1799,85 +1768,20 @@ namespace detail { return residualNorms; } - template - void - FullyImplicitBlackoilSolver::detectNewtonOscillations(const std::vector>& 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& F0 = residual_history[it]; - const std::vector& F1 = residual_history[it - 1]; - const std::vector& 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 - void - FullyImplicitBlackoilSolver::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 + template double - FullyImplicitBlackoilSolver::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 + BlackoilModelBase::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 @@ -1948,9 +1852,12 @@ namespace detail { } } - template + + + + template bool - FullyImplicitBlackoilSolver::getConvergence(const double dt, const int iteration) + BlackoilModelBase::getConvergence(const double dt, const int iteration) { const double tol_mb = param_.tolerance_mb_; const double tol_cnv = param_.tolerance_cnv_; @@ -2051,9 +1958,11 @@ namespace detail { } - template + + + template ADB - FullyImplicitBlackoilSolver::fluidViscosity(const int phase, + BlackoilModelBase::fluidViscosity(const int phase, const ADB& p , const ADB& temp , const ADB& rs , @@ -2078,9 +1987,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::fluidReciprocFVF(const int phase, + BlackoilModelBase::fluidReciprocFVF(const int phase, const ADB& p , const ADB& temp , const ADB& rs , @@ -2105,9 +2014,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::fluidDensity(const int phase, + BlackoilModelBase::fluidDensity(const int phase, const ADB& p , const ADB& temp , const ADB& rs , @@ -2133,9 +2042,9 @@ namespace detail { - template + template V - FullyImplicitBlackoilSolver::fluidRsSat(const V& p, + BlackoilModelBase::fluidRsSat(const V& p, const V& satOil, const std::vector& cells) const { @@ -2146,9 +2055,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::fluidRsSat(const ADB& p, + BlackoilModelBase::fluidRsSat(const ADB& p, const ADB& satOil, const std::vector& cells) const { @@ -2156,9 +2065,9 @@ namespace detail { } - template + template V - FullyImplicitBlackoilSolver::fluidRvSat(const V& p, + BlackoilModelBase::fluidRvSat(const V& p, const V& satOil, const std::vector& cells) const { @@ -2169,9 +2078,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::fluidRvSat(const ADB& p, + BlackoilModelBase::fluidRvSat(const ADB& p, const ADB& satOil, const std::vector& cells) const { @@ -2180,9 +2089,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::poroMult(const ADB& p) const + BlackoilModelBase::poroMult(const ADB& p) const { const int n = p.size(); if (rock_comp_props_ && rock_comp_props_->isActive()) { @@ -2208,9 +2117,9 @@ namespace detail { - template + template ADB - FullyImplicitBlackoilSolver::transMult(const ADB& p) const + BlackoilModelBase::transMult(const ADB& p) const { const int n = p.size(); if (rock_comp_props_ && rock_comp_props_->isActive()) { @@ -2234,9 +2143,9 @@ namespace detail { /* - template + template void - FullyImplicitBlackoilSolver:: + BlackoilModelBase:: classifyCondition(const SolutionState& state, std::vector& cond ) const { @@ -2276,9 +2185,9 @@ namespace detail { } */ - template + template void - FullyImplicitBlackoilSolver::classifyCondition(const BlackoilState& state) + BlackoilModelBase::classifyCondition(const BlackoilState& state) { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -2314,9 +2223,9 @@ namespace detail { } - template + template void - FullyImplicitBlackoilSolver::updatePrimalVariableFromState(const BlackoilState& state) + BlackoilModelBase::updatePrimalVariableFromState(const BlackoilState& state) { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -2364,9 +2273,9 @@ namespace detail { /// Update the phaseCondition_ member based on the primalVariable_ member. - template + template void - FullyImplicitBlackoilSolver::updatePhaseCondFromPrimalVariable() + BlackoilModelBase::updatePhaseCondFromPrimalVariable() { if (! active_[Gas]) { OPM_THROW(std::logic_error, "updatePhaseCondFromPrimarVariable() logic requires active gas phase."); @@ -2396,3 +2305,5 @@ namespace detail { } // namespace Opm + +#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED