From 26484e91a590d592974f699a2de0faf5fdf6753e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 24 May 2015 20:38:07 +0200 Subject: [PATCH] Transform BlackoilPolymerModel to inherit BlackoilModelBase. The class still contains surplus implementations though. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 272 ++++----- .../BlackoilPolymerModel_impl.hpp | 519 ++---------------- 2 files changed, 147 insertions(+), 644 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index f36ac964a..ff0afc342 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -21,31 +21,15 @@ #ifndef OPM_BLACKOILPOLYMERMODEL_HEADER_INCLUDED #define OPM_BLACKOILPOLYMERMODEL_HEADER_INCLUDED -#include - -#include -#include -#include -#include -#include +#include +#include #include #include - -#include - -struct UnstructuredGrid; -struct Wells; +#include +#include namespace Opm { - namespace parameter { class ParameterGroup; } - class DerivedGeology; - class RockCompressibility; - class NewtonIterationBlackoilInterface; - class PolymerBlackoilState; - class WellStateFullyImplicitBlackoilPolymer; - - /// A model implementation for three-phase black oil with polymer. /// /// The simulator is capable of handling three-phase problems @@ -56,36 +40,15 @@ namespace Opm { /// It uses automatic differentiation via the class AutoDiffBlock /// to simplify assembly of the jacobian matrix. template - class BlackoilPolymerModel + class BlackoilPolymerModel : public BlackoilModelBase > { public: // --------- Types and enums --------- - typedef AutoDiffBlock ADB; - typedef ADB::V V; - typedef ADB::M M; - typedef PolymerBlackoilState ReservoirState; - typedef WellStateFullyImplicitBlackoilPolymer 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 --------- + typedef BlackoilModelBase > Base; + typedef typename Base::ReservoirState ReservoirState; + typedef typename Base::WellState WellState; /// Construct the model. It will retain references to the /// arguments of this functions, and they are expected to @@ -101,28 +64,18 @@ namespace Opm { /// \param[in] has_vapoil turn on vaporized oil feature /// \param[in] has_polymer turn on polymer feature /// \param[in] terminal_output request output to cout/cerr - BlackoilPolymerModel(const ModelParameters& param, - const Grid& grid , - const BlackoilPropsAdInterface& fluid, - const DerivedGeology& geo , - const RockCompressibility* rock_comp_props, - const PolymerPropsAd& polymer_props_ad, - const Wells* wells, + BlackoilPolymerModel(const typename Base::ModelParameters& param, + const Grid& grid, + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo, + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells* wells, const NewtonIterationBlackoilInterface& linsolver, - const bool has_disgas, - const bool has_vapoil, - const bool has_polymer, - const bool terminal_output); - - /// \brief Set threshold pressures that prevent or reduce flow. - /// This prevents flow across faces if the potential - /// difference is less than the threshold. If the potential - /// difference is greater, the threshold value is subtracted - /// before calculating flow. This is treated symmetrically, so - /// flow is prevented or reduced in both directions equally. - /// \param[in] threshold_pressures_by_face array of size equal to the number of faces - /// of the grid passed in the constructor. - void setThresholdPressures(const std::vector& threshold_pressures_by_face); + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + const bool terminal_output); /// Called once before each time step. /// \param[in] dt time step size @@ -147,19 +100,15 @@ namespace Opm { void assemble(const ReservoirState& reservoir_state, WellState& well_state, const bool initial_assembly); - + // void assemble(const PolymerBlackoilState& reservoir_state, + // WellStateFullyImplicitBlackoilPolymer& 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; @@ -172,138 +121,86 @@ namespace Opm { 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; - - private: + protected: // --------- Types and enums --------- - typedef Eigen::Array DataBlock; - - struct ReservoirResidualQuant { - ReservoirResidualQuant(); - std::vector accum; // Accumulations - ADB mflux; // Mass flux (surface conditions) - ADB b; // Reciprocal FVF - ADB head; // Pressure drop across int. interfaces - ADB mob; // Phase mobility (per cell) - }; - - struct SolutionState { - SolutionState(const int np); - ADB pressure; - ADB temperature; - std::vector saturation; - ADB rs; - ADB rv; - ADB concentration; - ADB qs; - ADB bhp; - // Below are quantities stored in the state for optimization purposes. - std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active. - }; - - struct WellOps { - WellOps(const Wells* wells); - M w2p; // well -> perf (scatter) - M p2w; // perf -> well (gather) - }; - - enum { Water = BlackoilPropsAdInterface::Water, - Oil = BlackoilPropsAdInterface::Oil , - Gas = BlackoilPropsAdInterface::Gas , - MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases - }; - - enum PrimalVariables { Sg = 0, RS = 1, RV = 2 }; + typedef typename Base::SolutionState SolutionState; + typedef typename Base::DataBlock DataBlock; // --------- Data members --------- - const Grid& grid_; - const BlackoilPropsAdInterface& fluid_; - const DerivedGeology& geo_; - const RockCompressibility* rock_comp_props_; - const PolymerPropsAd& polymer_props_ad_; - const Wells* wells_; - const NewtonIterationBlackoilInterface& linsolver_; - // For each canonical phase -> true if active - const std::vector active_; - // Size = # active phases. Maps active -> canonical phase indices. - const std::vector canph_; - const std::vector cells_; // All grid cells - HelperOps ops_; - const WellOps wops_; - V cmax_; - const bool has_disgas_; - const bool has_vapoil_; + const PolymerPropsAd& polymer_props_ad_; const bool has_polymer_; const int poly_pos_; + V cmax_; - ModelParameters param_; - bool use_threshold_pressure_; - V threshold_pressures_by_interior_face_; + // Need to declare Base members we want to use here. + using Base::grid_; + using Base::fluid_; + using Base::geo_; + using Base::rock_comp_props_; + using Base::wells_; + using Base::linsolver_; + using Base::active_; + using Base::canph_; + using Base::cells_; + using Base::ops_; + using Base::wops_; + using Base::has_disgas_; + using Base::has_vapoil_; + using Base::param_; + using Base::use_threshold_pressure_; + using Base::threshold_pressures_by_interior_face_; + using Base::rq_; + using Base::phaseCondition_; + using Base::well_perforation_pressure_diffs_; + using Base::residual_; + using Base::terminal_output_; + using Base::primalVariable_; + using Base::pvdt_; - std::vector rq_; - std::vector phaseCondition_; - V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. + // --------- Protected methods --------- - LinearisedBlackoilResidual residual_; - - /// \brief Whether we print something to std::cout - bool terminal_output_; - - std::vector primalVariable_; - V pvdt_; - - // --------- Private methods --------- - - // return true if wells are available - bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; } - // return wells object - const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; } + // Need to declare Base members we want to use here. + using Base::wellsActive; + using Base::wells; SolutionState - constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoilPolymer& xw) const; + constantState(const ReservoirState& x, + const WellState& xw) const; void makeConstantState(SolutionState& state) const; SolutionState - variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoilPolymer& xw) const; + variableState(const ReservoirState& x, + const WellState& xw) const; void computeAccum(const SolutionState& state, const int aix ); void computeWellConnectionPressures(const SolutionState& state, - const WellStateFullyImplicitBlackoilPolymer& xw); + const WellState& xw); void addWellControlEq(const SolutionState& state, - const WellStateFullyImplicitBlackoilPolymer& xw, + const WellState& xw, const V& aliveWells); void addWellEq(const SolutionState& state, - WellStateFullyImplicitBlackoilPolymer& xw, + WellState& xw, V& aliveWells); - void updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const; + void updateWellControls(WellState& xw) const; std::vector computePressures(const SolutionState& state) const; @@ -331,7 +228,7 @@ namespace Opm { const SolutionState& state ); void - computeCmax(PolymerBlackoilState& state); + computeCmax(ReservoirState& state); ADB computeMc(const SolutionState& state) const; @@ -396,16 +293,16 @@ namespace Opm { std::vector& cond ) const; const std::vector - phaseCondition() const {return phaseCondition_;} + phaseCondition() const {return this->phaseCondition_;} void - classifyCondition(const PolymerBlackoilState& state); + classifyCondition(const ReservoirState& state); /// update the primal variable for Sg, Rv or Rs. The Gas phase must /// be active to call this method. void - updatePrimalVariableFromState(const PolymerBlackoilState& state); + updatePrimalVariableFromState(const ReservoirState& state); /// Update the phaseCondition_ member based on the primalVariable_ member. void @@ -441,12 +338,39 @@ namespace Opm { 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_; } + double dpMaxRel() const { return this->param_.dp_max_rel_; } + double dsMax() const { return this->param_.ds_max_; } + double drMaxRel() const { return this->param_.dr_max_rel_; } + double maxResidualAllowed() const { return this->param_.max_residual_allowed_; } }; + + + + /// Need to include concentration in our state variables, otherwise all is as + /// the default blackoil model. + struct BlackoilPolymerSolutionState : public DefaultBlackoilSolutionState + { + explicit BlackoilPolymerSolutionState(const int np) + : DefaultBlackoilSolutionState(np), + concentration( ADB::null()) + { + } + ADB concentration; + }; + + + + /// Providing types by template specialisation of ModelTraits for BlackoilPolymerModel. + template + struct ModelTraits< BlackoilPolymerModel > + { + typedef PolymerBlackoilState ReservoirState; + typedef WellStateFullyImplicitBlackoilPolymer WellState; + typedef BlackoilModelParameters ModelParameters; + typedef BlackoilPolymerSolutionState SolutionState; + }; + } // namespace Opm #include "BlackoilPolymerModel_impl.hpp" diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c4b1cba23..c0285914b 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -25,8 +25,6 @@ #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED #include -#include -#include #include #include @@ -50,205 +48,57 @@ #include #include #include -//#include - -// A debugging utility. -#define OPM_AD_DUMP(foo) \ - do { \ - std::cout << "==========================================\n" \ - << #foo ":\n" \ - << collapseJacs(foo) << std::endl; \ - } while (0) - -#define OPM_AD_DUMPVAL(foo) \ - do { \ - std::cout << "==========================================\n" \ - << #foo ":\n" \ - << foo.value() << std::endl; \ - } while (0) - -#define OPM_AD_DISKVAL(foo) \ - do { \ - std::ofstream os(#foo); \ - os.precision(16); \ - os << foo.value() << std::endl; \ - } while (0) - namespace Opm { -typedef AutoDiffBlock ADB; -typedef ADB::V V; -typedef ADB::M M; -typedef Eigen::Array DataBlock; -namespace detail { + namespace detail { - - std::vector - buildAllCells(const int nc) - { - std::vector all_cells(nc); - - for (int c = 0; c < nc; ++c) { all_cells[c] = c; } - - return all_cells; - } - - - - template - std::vector - activePhases(const PU& pu) - { - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - std::vector active(maxnp, false); - - for (int p = 0; p < pu.MaxNumPhases; ++p) { - active[ p ] = pu.phase_used[ p ] != 0; - } - - return active; - } - - - - template - std::vector - active2Canonical(const PU& pu) - { - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - std::vector act2can(maxnp, -1); - - for (int phase = 0; phase < maxnp; ++phase) { - if (pu.phase_used[ phase ]) { - act2can[ pu.phase_pos[ phase ] ] = phase; + template + int polymerPos(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + int pos = 0; + for (int phase = 0; phase < maxnp; ++phase) { + if (pu.phase_used[phase]) { + pos++; + } } + + return pos; } - return act2can; - } + } // namespace detail - template - int polymerPos(const PU& pu) - { - const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - int pos = 0; - for (int phase = 0; phase < maxnp; ++phase) { - if (pu.phase_used[phase]) { - pos++; - } - } - - return pos; - } - -} // namespace detail - - template - void BlackoilPolymerModel::ModelParameters:: - reset() - { - // default values for the solver parameters - dp_max_rel_ = 1.0e9; - ds_max_ = 0.2; - dr_max_rel_ = 1.0e9; - max_residual_allowed_ = 1e7; - tolerance_mb_ = 1.0e-5; - tolerance_cnv_ = 1.0e-2; - tolerance_wells_ = 5.0e-1; - } - - template - BlackoilPolymerModel::ModelParameters:: - ModelParameters() - { - // set default values - reset(); - } - - template - BlackoilPolymerModel::ModelParameters:: - ModelParameters( const parameter::ParameterGroup& param ) - { - // set default values - reset(); - - // overload with given parameters - dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_); - ds_max_ = param.getDefault("ds_max", ds_max_); - dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_); - max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_); - tolerance_mb_ = param.getDefault("tolerance_mb", tolerance_mb_); - tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_); - tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ ); - } - - template - BlackoilPolymerModel:: - BlackoilPolymerModel(const ModelParameters& param, - const Grid& grid , - const BlackoilPropsAdInterface& fluid, - const DerivedGeology& geo , - const RockCompressibility* rock_comp_props, - const PolymerPropsAd& polymer_props_ad, - const Wells* wells, - const NewtonIterationBlackoilInterface& linsolver, - const bool has_disgas, - const bool has_vapoil, - const bool has_polymer, - const bool terminal_output) - : grid_ (grid) - , fluid_ (fluid) - , geo_ (geo) - , rock_comp_props_(rock_comp_props) - , polymer_props_ad_(polymer_props_ad) - , wells_ (wells) - , linsolver_ (linsolver) - , active_(detail::activePhases(fluid.phaseUsage())) - , canph_ (detail::active2Canonical(fluid.phaseUsage())) - , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid))) - , ops_ (grid) - , wops_ (wells_) - , cmax_(V::Zero(Opm::AutoDiffGrid::numCells(grid))) - , has_disgas_(has_disgas) - , has_vapoil_(has_vapoil) - , has_polymer_(has_polymer) - , poly_pos_(detail::polymerPos(fluid.phaseUsage())) - , param_( param ) - , use_threshold_pressure_(false) - , rq_ (fluid.numPhases()) - , phaseCondition_(AutoDiffGrid::numCells(grid)) - , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), - ADB::null(), - ADB::null() } ) - , terminal_output_ (terminal_output) + BlackoilPolymerModel::BlackoilPolymerModel(const typename Base::ModelParameters& param, + const Grid& grid, + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo, + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells* wells, + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + const bool terminal_output) + : Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver, + has_disgas, has_vapoil, terminal_output), + polymer_props_ad_(polymer_props_ad), + has_polymer_(has_polymer), + poly_pos_(detail::polymerPos(fluid.phaseUsage())) { -#if HAVE_MPI - if ( terminal_output_ ) { - if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& info = - boost::any_cast(linsolver_.parallelInformation()); - // Only rank 0 does print to std::cout if terminal_output is enabled - terminal_output_ = (info.communicator().rank()==0); - } - } -#endif if (has_polymer_) { if (!active_[Water]) { OPM_THROW(std::logic_error, "Polymer must solved in water!\n"); } // If deck has polymer, residual_ should contain polymer equation. - rq_.resize(fluid_.numPhases()+1); - residual_.material_balance_eq.resize(fluid_.numPhases()+1, ADB::null()); + rq_.resize(fluid_.numPhases() + 1); + residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null()); assert(poly_pos_ == fluid_.numPhases()); } } @@ -260,12 +110,9 @@ namespace detail { BlackoilPolymerModel:: prepareStep(const double dt, ReservoirState& reservoir_state, - WellState& /* well_state */) + WellState& well_state) { - pvdt_ = geo_.poreVolume() / dt; - if (active_[Gas]) { - updatePrimalVariableFromState(reservoir_state); - } + Base::prepareStep(dt, reservoir_state, well_state); // Initial max concentration of this time step from PolymerBlackoilState. cmax_ = Eigen::Map(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_)); } @@ -286,144 +133,11 @@ namespace detail { - template - int - BlackoilPolymerModel:: - sizeNonLinear() const - { - return residual_.sizeNonLinear(); - } - - - - - template - int - BlackoilPolymerModel:: - linearIterationsLastSolve() const - { - return linsolver_.iterations(); - } - - - - - template - bool - BlackoilPolymerModel:: - terminalOutputEnabled() const - { - return terminal_output_; - } - - - - - template - int - BlackoilPolymerModel:: - numPhases() const - { - return fluid_.numPhases(); - } - - - - - template - void - BlackoilPolymerModel:: - setThresholdPressures(const std::vector& threshold_pressures_by_face) - { - const int num_faces = AutoDiffGrid::numFaces(grid_); - if (int(threshold_pressures_by_face.size()) != num_faces) { - OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces."); - } - use_threshold_pressure_ = true; - // Map to interior faces. - const int num_ifaces = ops_.internal_faces.size(); - threshold_pressures_by_interior_face_.resize(num_ifaces); - for (int ii = 0; ii < num_ifaces; ++ii) { - threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]]; - } - } - - - - - - template - BlackoilPolymerModel::ReservoirResidualQuant::ReservoirResidualQuant() - : accum(2, ADB::null()) - , mflux( ADB::null()) - , b ( ADB::null()) - , head ( ADB::null()) - , mob ( ADB::null()) - { - } - - - - - - template - BlackoilPolymerModel::SolutionState::SolutionState(const int np) - : pressure ( ADB::null()) - , temperature( ADB::null()) - , saturation(np, ADB::null()) - , rs ( ADB::null()) - , rv ( ADB::null()) - , concentration( ADB::null()) - , qs ( ADB::null()) - , bhp ( ADB::null()) - , canonical_phase_pressures(3, ADB::null()) - { - } - - - - - - template - BlackoilPolymerModel:: - WellOps::WellOps(const Wells* wells) - : w2p(), - p2w() - { - if( wells ) - { - w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); - p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); - - const int nw = wells->number_of_wells; - const int* const wpos = wells->well_connpos; - - typedef Eigen::Triplet Tri; - - std::vector scatter, gather; - scatter.reserve(wpos[nw]); - gather .reserve(wpos[nw]); - - for (int w = 0, i = 0; w < nw; ++w) { - for (; i < wpos[ w + 1 ]; ++i) { - scatter.push_back(Tri(i, w, 1.0)); - gather .push_back(Tri(w, i, 1.0)); - } - } - - w2p.setFromTriplets(scatter.begin(), scatter.end()); - p2w.setFromTriplets(gather .begin(), gather .end()); - } - } - - - - template typename BlackoilPolymerModel::SolutionState - BlackoilPolymerModel::constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoilPolymer& xw) const + BlackoilPolymerModel::constantState(const ReservoirState& x, + const WellState& xw) const { auto state = variableState(x, xw); makeConstantState(state); @@ -466,8 +180,8 @@ namespace detail { template typename BlackoilPolymerModel::SolutionState - BlackoilPolymerModel::variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoilPolymer& xw) const + BlackoilPolymerModel::variableState(const ReservoirState& x, + const WellState& xw) const { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -693,7 +407,7 @@ namespace detail { template - void BlackoilPolymerModel::computeCmax(PolymerBlackoilState& state) + void BlackoilPolymerModel::computeCmax(ReservoirState& state) { const int nc = AutoDiffGrid::numCells(grid_); V tmp = V::Zero(nc); @@ -709,7 +423,7 @@ namespace detail { template void BlackoilPolymerModel::computeWellConnectionPressures(const SolutionState& state, - const WellStateFullyImplicitBlackoilPolymer& xw) + const WellState& xw) { if( ! wellsActive() ) return ; @@ -803,8 +517,8 @@ namespace detail { template void BlackoilPolymerModel:: - assemble(const PolymerBlackoilState& reservoir_state, - WellStateFullyImplicitBlackoilPolymer& well_state, + assemble(const ReservoirState& reservoir_state, + WellState& well_state, const bool initial_assembly) { using namespace Opm::AutoDiffGrid; @@ -898,7 +612,7 @@ namespace detail { template void BlackoilPolymerModel::addWellEq(const SolutionState& state, - WellStateFullyImplicitBlackoilPolymer& xw, + WellState& xw, V& aliveWells) { if( ! wellsActive() ) return ; @@ -1073,85 +787,8 @@ namespace detail { - namespace detail - { - double rateToCompare(const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const double* distr) - { - double rate = 0.0; - for (int phase = 0; phase < num_phases; ++phase) { - // Important: well_phase_flow_rate is ordered with all phase rates for first - // well first, then all phase rates for second well etc. - rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase]; - } - return rate; - } - - bool constraintBroken(const std::vector& bhp, - const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const WellType& well_type, - const WellControls* wc, - const int ctrl_index) - { - const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); - const double target = well_controls_iget_target(wc, ctrl_index); - const double* distr = well_controls_iget_distr(wc, ctrl_index); - - bool broken = false; - - switch (well_type) { - case INJECTOR: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] > target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) > target; - break; - } - } - break; - - case PRODUCER: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] < target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - // Note that the rates compared below are negative, - // so breaking the constraints means: too high flow rate - // (as for injection). - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) < target; - break; - } - } - break; - - default: - OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); - } - - return broken; - } - } // namespace detail - - - - template - void BlackoilPolymerModel::updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const + void BlackoilPolymerModel::updateWellControls(WellState& xw) const { if( ! wellsActive() ) return ; @@ -1227,7 +864,7 @@ namespace detail { template void BlackoilPolymerModel::addWellControlEq(const SolutionState& state, - const WellStateFullyImplicitBlackoilPolymer& xw, + const WellState& xw, const V& aliveWells) { if( ! wellsActive() ) return; @@ -1306,68 +943,10 @@ namespace detail { - namespace detail - { - /// \brief Compute the L-infinity norm of a vector - /// \warn This function is not suitable to compute on the well equations. - /// \param a The container to compute the infinity norm on. - /// It has to have one entry for each cell. - /// \param info In a parallel this holds the information about the data distribution. - double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() ) - { -#if HAVE_MPI - if ( pinfo.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& real_info = - boost::any_cast(pinfo); - double result=0; - real_info.computeReduction(a.value(), Reduction::makeGlobalMaxFunctor(), result); - return result; - } - else -#endif - { - if( a.value().size() > 0 ) { - return a.value().matrix().lpNorm (); - } - else { // this situation can occur when no wells are present - return 0.0; - } - } - (void)pinfo; // Suppress unused warning for non-MPI. - } - - /// \brief Compute the L-infinity norm of a vector representing a well equation. - /// \param a The container to compute the infinity norm on. - /// \param info In a parallel this holds the information about the data distribution. - double infinityNormWell( const ADB& a, const boost::any& pinfo ) - { - double result=0; - if( a.value().size() > 0 ) { - result = a.value().matrix().lpNorm (); - } -#if HAVE_MPI - if ( pinfo.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& real_info = - boost::any_cast(pinfo); - result = real_info.communicator().max(result); - } -#endif - return result; - (void)pinfo; // Suppress unused warning for non-MPI. - } - - } // namespace detail - - - - - template void BlackoilPolymerModel::updateState(const V& dx, - PolymerBlackoilState& reservoir_state, - WellStateFullyImplicitBlackoilPolymer& well_state) + ReservoirState& reservoir_state, + WellState& well_state) { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); @@ -2321,7 +1900,7 @@ namespace detail { template void - BlackoilPolymerModel::classifyCondition(const PolymerBlackoilState& state) + BlackoilPolymerModel::classifyCondition(const ReservoirState& state) { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -2359,7 +1938,7 @@ namespace detail { template void - BlackoilPolymerModel::updatePrimalVariableFromState(const PolymerBlackoilState& state) + BlackoilPolymerModel::updatePrimalVariableFromState(const ReservoirState& state) { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_);