diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index c3b502c41..8b616c0af 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -43,15 +43,15 @@ namespace Opm { class RockCompressibility; class NewtonIterationBlackoilInterface; class PolymerBlackoilState; - class WellStateFullyImplicitBlackoil; + class WellStateFullyImplicitBlackoilPolymer; - /// A fully implicit solver for the black-oil-polymer problem. + /// A model implementation for three-phase black oil with polymer. /// /// The simulator is capable of handling three-phase problems - /// where gas can be dissolved in oil (but not vice versa). It - /// uses an industry-standard TPFA discretization with per-phase - /// upwind weighting of mobilities. + /// where gas can be dissolved in oil and vice versa, with polymer + /// in the water phase. 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. @@ -59,33 +59,35 @@ namespace Opm { class BlackoilPolymerModel { public: - // the Newton relaxation type - enum RelaxType { DAMPEN, SOR }; - // class holding the solver parameters - struct SolverParameter + // --------- 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 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 + 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(); + ModelParameter( const parameter::ParameterGroup& param ); + ModelParameter(); void reset(); }; - /// 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 @@ -95,7 +97,11 @@ namespace Opm { /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] wells well structure /// \param[in] linsolver linear solver - BlackoilPolymerModel(const SolverParameter& param, + /// \param[in] has_disgas turn on dissolved gas + /// \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 ModelParameter& param, const Grid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , @@ -118,30 +124,70 @@ 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 , - PolymerBlackoilState& state , - WellStateFullyImplicitBlackoil& wstate, - const std::vector& polymer_inflow); + /// 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, + ReservoirState& reservoir_state, + WellState& well_state); - unsigned int newtonIterations () const { return newtonIterations_; } - unsigned int linearIterations () const { return linearIterations_; } + /// Called once after each time step. + /// \param[in] dt time step size + /// \param[in] reservoir_state reservoir state variables + /// \param[in] well_state well state variables + void afterStep(const double dt, + ReservoirState& reservoir_state, + WellState& well_state); + + /// Assemble the residual and Jacobian of the nonlinear system. + /// \param[in] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep + void assemble(const ReservoirState& reservoir_state, + WellState& well_state, + const bool initial_assembly); + + /// \brief Compute the residual norms of the mass balance for each phase, + /// the well flux, and the well equation. + /// \return a vector that contains for each phase the norm of the mass balance + /// and afterwards the norm of the residual of the well flux and the well equation. + std::vector computeResidualNorms() const; + + /// The size (number of unknowns) of the nonlinear system of equations. + int sizeNonLinear() const; + + /// Number of linear iterations used in last call to solveJacobianSystem(). + int linearIterationsLastSolve() const; + + /// Solve the Jacobian system Jx = r where J is the Jacobian and + /// r is the residual. + V solveJacobianSystem() const; + + /// Apply an update to the primary variables, chopped if appropriate. + /// \param[in] dx updates to apply to primary variables + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + void updateState(const V& dx, + ReservoirState& reservoir_state, + WellState& well_state); + + /// Return true if output to cout is wanted. + bool terminalOutput() 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 ; } @@ -231,47 +279,33 @@ namespace Opm { SolutionState constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) const; + const WellStateFullyImplicitBlackoilPolymer& xw) const; void makeConstantState(SolutionState& state) const; SolutionState variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) const; + const WellStateFullyImplicitBlackoilPolymer& xw) const; void computeAccum(const SolutionState& state, const int aix ); void computeWellConnectionPressures(const SolutionState& state, - const WellStateFullyImplicitBlackoil& xw); + const WellStateFullyImplicitBlackoilPolymer& xw); void addWellControlEq(const SolutionState& state, - const WellStateFullyImplicitBlackoil& xw, + const WellStateFullyImplicitBlackoilPolymer& xw, const V& aliveWells); void addWellEq(const SolutionState& state, - WellStateFullyImplicitBlackoil& xw, - V& aliveWells, - const std::vector& polymer_inflow); + WellStateFullyImplicitBlackoilPolymer& xw, + V& aliveWells); - void updateWellControls(WellStateFullyImplicitBlackoil& xw) const; - - void - assemble(const V& dtpv, - const PolymerBlackoilState& x, - const bool initial_assembly, - WellStateFullyImplicitBlackoil& xw, - const std::vector& polymer_inflow); - - V solveJacobianSystem() const; - - void updateState(const V& dx, - PolymerBlackoilState& state, - WellStateFullyImplicitBlackoil& well_state); + void updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const; std::vector computePressures(const SolutionState& state) const; @@ -306,15 +340,6 @@ namespace Opm { void applyThresholdPressures(ADB& dp); - double - residualNorm() const; - - /// \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 , @@ -388,10 +413,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 @@ -420,21 +441,9 @@ namespace Opm { std::array& B_avg, int nc) 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_; } }; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 96be869d9..b05a97aea 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -32,7 +33,6 @@ #include #include #include -#include #include #include @@ -151,19 +151,13 @@ namespace detail { } // namespace detail template - void BlackoilPolymerModel::SolverParameter:: + void BlackoilPolymerModel::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; @@ -172,16 +166,16 @@ namespace detail { } template - BlackoilPolymerModel::SolverParameter:: - SolverParameter() + BlackoilPolymerModel::ModelParameter:: + ModelParameter() { // set default values reset(); } template - BlackoilPolymerModel::SolverParameter:: - SolverParameter( const parameter::ParameterGroup& param ) + BlackoilPolymerModel::ModelParameter:: + ModelParameter( const parameter::ParameterGroup& param ) { // set default values reset(); @@ -190,40 +184,27 @@ 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 BlackoilPolymerModel:: - BlackoilPolymerModel(const SolverParameter& 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) + BlackoilPolymerModel(const ModelParameter& 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) @@ -276,6 +257,81 @@ namespace detail { + template + void + BlackoilPolymerModel:: + prepareStep(const double dt, + ReservoirState& reservoir_state, + WellState& /* well_state */) + { + pvdt_ = geo_.poreVolume() / dt; + if (active_[Gas]) { + updatePrimalVariableFromState(reservoir_state); + } + // Initial max concentration of this time step from PolymerBlackoilState. + cmax_ = Eigen::Map(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_)); + } + + + + + template + void + BlackoilPolymerModel:: + afterStep(const double /* dt */, + ReservoirState& reservoir_state, + WellState& /* well_state */) + { + computeCmax(reservoir_state); + } + + + + + template + int + BlackoilPolymerModel:: + sizeNonLinear() const + { + return residual_.sizeNonLinear(); + } + + + + + template + int + BlackoilPolymerModel:: + linearIterationsLastSolve() const + { + return linsolver_.iterations(); + } + + + + + template + bool + BlackoilPolymerModel:: + terminalOutput() const + { + return terminal_output_; + } + + + + + template + int + BlackoilPolymerModel:: + numPhases() const + { + return fluid_.numPhases(); + } + + + + template void BlackoilPolymerModel:: @@ -297,89 +353,6 @@ namespace detail { - template - int - BlackoilPolymerModel:: - step(const double dt, - PolymerBlackoilState& x , - WellStateFullyImplicitBlackoil& xw, - const std::vector& polymer_inflow) - { - const V pvdt = geo_.poreVolume() / dt; - - // Initial max concentration of this time step from PolymerBlackoilState. - cmax_ = Eigen::Map(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_)); - 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, polymer_inflow); - - - 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, polymer_inflow); - - 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; - // Update max concentration. - computeCmax(x); - return linearIterations; - } - - - - template BlackoilPolymerModel::ReservoirResidualQuant::ReservoirResidualQuant() @@ -452,7 +425,7 @@ namespace detail { template typename BlackoilPolymerModel::SolutionState BlackoilPolymerModel::constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) const + const WellStateFullyImplicitBlackoilPolymer& xw) const { auto state = variableState(x, xw); makeConstantState(state); @@ -496,7 +469,7 @@ namespace detail { template typename BlackoilPolymerModel::SolutionState BlackoilPolymerModel::variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) const + const WellStateFullyImplicitBlackoilPolymer& xw) const { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); @@ -738,7 +711,7 @@ namespace detail { template void BlackoilPolymerModel::computeWellConnectionPressures(const SolutionState& state, - const WellStateFullyImplicitBlackoil& xw) + const WellStateFullyImplicitBlackoilPolymer& xw) { if( ! wellsActive() ) return ; @@ -832,11 +805,9 @@ namespace detail { template void BlackoilPolymerModel:: - assemble(const V& pvdt, - const PolymerBlackoilState& x , - const bool initial_assembly, - WellStateFullyImplicitBlackoil& xw, - const std::vector& polymer_inflow) + assemble(const PolymerBlackoilState& x, + WellStateFullyImplicitBlackoilPolymer& xw, + const bool initial_assembly) { using namespace Opm::AutoDiffGrid; @@ -883,7 +854,7 @@ namespace detail { for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); residual_.material_balance_eq[ phaseIdx ] = - pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + ops_.div*rq_[phaseIdx].mflux; } @@ -913,13 +884,13 @@ namespace detail { // Add polymer equation. if (has_polymer_) { - residual_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0]) + residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0]) + ops_.div*rq_[poly_pos_].mflux; } // Add contribution from wells and set up the well equations. V aliveWells; - addWellEq(state, xw, aliveWells, polymer_inflow); + addWellEq(state, xw, aliveWells); addWellControlEq(state, xw, aliveWells); } @@ -929,9 +900,8 @@ namespace detail { template void BlackoilPolymerModel::addWellEq(const SolutionState& state, - WellStateFullyImplicitBlackoil& xw, - V& aliveWells, - const std::vector& polymer_inflow) + WellStateFullyImplicitBlackoilPolymer& xw, + V& aliveWells) { if( ! wellsActive() ) return ; @@ -1065,7 +1035,7 @@ namespace detail { // Add well contributions to polymer mass balance equation if (has_polymer_) { const ADB mc = computeMc(state); - const V polyin = Eigen::Map(&polymer_inflow[0], nc); + const V polyin = Eigen::Map(xw.polymerInflow().data(), nc); const V poly_in_perf = subset(polyin, well_cells); const V poly_mc_perf = subset(mc, well_cells).value(); residual_.material_balance_eq[poly_pos_] -= superset(cq_ps[pu.phase_pos[Water]] * poly_mc_perf @@ -1183,7 +1153,7 @@ namespace detail { template - void BlackoilPolymerModel::updateWellControls(WellStateFullyImplicitBlackoil& xw) const + void BlackoilPolymerModel::updateWellControls(WellStateFullyImplicitBlackoilPolymer& xw) const { if( ! wellsActive() ) return ; @@ -1259,8 +1229,8 @@ namespace detail { template void BlackoilPolymerModel::addWellControlEq(const SolutionState& state, - const WellStateFullyImplicitBlackoil& xw, - const V& aliveWells) + const WellStateFullyImplicitBlackoilPolymer& xw, + const V& aliveWells) { if( ! wellsActive() ) return; @@ -1359,8 +1329,8 @@ namespace detail { template void BlackoilPolymerModel::updateState(const V& dx, - PolymerBlackoilState& state, - WellStateFullyImplicitBlackoil& well_state) + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoilPolymer& well_state) { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); @@ -1835,30 +1805,6 @@ namespace detail { - template - double - BlackoilPolymerModel::residualNorm() const - { - double globalNorm = 0; - std::vector::const_iterator quantityIt = residual_.material_balance_eq.begin(); - const std::vector::const_iterator endQuantityIt = residual_.material_balance_eq.end(); - for (; quantityIt != endQuantityIt; ++quantityIt) { - const double quantityResid = (*quantityIt).value().matrix().norm(); - if (!std::isfinite(quantityResid)) { - const int trouble_phase = quantityIt - residual_.material_balance_eq.begin(); - OPM_THROW(Opm::NumericalProblem, - "Encountered a non-finite residual in material balance equation " - << trouble_phase); - } - globalNorm = std::max(globalNorm, quantityResid); - } - globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm()); - globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm()); - - return globalNorm; - } - - template std::vector BlackoilPolymerModel::computeResidualNorms() const @@ -1895,73 +1841,8 @@ namespace detail { return residualNorms; } - template - void - BlackoilPolymerModel::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 - BlackoilPolymerModel::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 double diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 53fc55672..73f987f49 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -30,7 +31,6 @@ #include #include -#include #include #include @@ -135,7 +135,7 @@ namespace Opm computeRESV(const std::size_t step, const Wells* wells, const BlackoilState& x, - WellStateFullyImplicitBlackoil& xw); + WellStateFullyImplicitBlackoilPolymer& xw); }; @@ -236,7 +236,7 @@ namespace Opm SimulatorReport SimulatorFullyImplicitBlackoilPolymer::Impl::run(SimulatorTimer& timer, PolymerBlackoilState& state) { - WellStateFullyImplicitBlackoil prev_well_state; + WellStateFullyImplicitBlackoilPolymer prev_well_state; // Create timers and file for writing timing info. Opm::time::StopWatch solver_timer; @@ -249,11 +249,10 @@ namespace Opm typedef T Grid; typedef BlackoilPolymerModel Model; - // typedef typename Model::ModelParameter ModelParam; - // ModelParam modelParam( param_ ); - // typedef NewtonSolver Solver; - // typedef typename Solver::SolverParameter SolverParam; - typedef typename Model::SolverParameter SolverParam; + typedef typename Model::ModelParameter ModelParam; + ModelParam modelParam( param_ ); + typedef NewtonSolver Solver; + typedef typename Solver::SolverParameter SolverParam; SolverParam solverParam( param_ ); //adaptive time stepping @@ -297,7 +296,7 @@ namespace Opm Opm::UgGridHelpers::beginFaceCentroids(grid_), props_.permeability()); const Wells* wells = wells_manager.c_wells(); - WellStateFullyImplicitBlackoil well_state; + WellStateFullyImplicitBlackoilPolymer well_state; well_state.init(wells, state.blackoilState(), prev_well_state); // compute polymer inflow @@ -316,7 +315,8 @@ namespace Opm polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), timer.simulationTimeElapsed() + timer.currentStepLength(), polymer_inflow_c); - + well_state.polymerInflow() = polymer_inflow_c; + // write simulation state at the report stage output_writer_.writeTimeStep( timer, state.blackoilState(), well_state ); @@ -330,10 +330,11 @@ namespace Opm // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - Model solver(solverParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_); + Model model(modelParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_); if (!threshold_pressures_by_face_.empty()) { - solver.setThresholdPressures(threshold_pressures_by_face_); + model.setThresholdPressures(threshold_pressures_by_face_); } + 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 @@ -344,7 +345,7 @@ namespace Opm // adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ ); // } else { // solve for complete report step - solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + solver.step(timer.currentStepLength(), state, well_state); // } // take time that was used to solve system for this reportStep @@ -511,7 +512,7 @@ namespace Opm Impl::computeRESV(const std::size_t step, const Wells* wells, const BlackoilState& x, - WellStateFullyImplicitBlackoil& xw) + WellStateFullyImplicitBlackoilPolymer& xw) { typedef SimFIBODetails::WellMap WellMap; diff --git a/opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp new file mode 100644 index 000000000..77327f15d --- /dev/null +++ b/opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp @@ -0,0 +1,39 @@ +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED +#define OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED + +#include + +namespace Opm +{ + + class WellStateFullyImplicitBlackoilPolymer : public WellStateFullyImplicitBlackoil + { + public: + std::vector& polymerInflow() { return polymer_inflow_; } + const std::vector& polymerInflow() const { return polymer_inflow_; } + private: + std::vector polymer_inflow_; + }; + +} // namespace Opm + +#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED