diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 9e28b85f1..0e27777fa 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -52,6 +52,10 @@ NEW_PROP_TAG(TolerancePressureMsWells); NEW_PROP_TAG(MaxPressureChangeMsWells); NEW_PROP_TAG(UseInnerIterationsMsWells); NEW_PROP_TAG(MaxInnerIterMsWells); +NEW_PROP_TAG(StrictInnerIterMsWells); +NEW_PROP_TAG(RelaxedFlowTolInnerIterMsw); +NEW_PROP_TAG(RelaxedPressureTolInnerIterMsw); +NEW_PROP_TAG(RegularizationFactorMsw); SET_SCALAR_PROP(FlowModelParameters, DbhpMaxRel, 1.0); SET_SCALAR_PROP(FlowModelParameters, DwellFractionMax, 0.2); @@ -73,8 +77,13 @@ SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01*1e5); SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 10*1e5); SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true); SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 100); +SET_INT_PROP(FlowModelParameters, StrictInnerIterMsWells, 40); +SET_SCALAR_PROP(FlowModelParameters, RegularizationFactorMsw, 1); SET_BOOL_PROP(FlowModelParameters, EnableWellOperabilityCheck, true); +SET_SCALAR_PROP(FlowModelParameters, RelaxedFlowTolInnerIterMsw, 1); +SET_SCALAR_PROP(FlowModelParameters, RelaxedPressureTolInnerIterMsw, 0.5e5); + // if openMP is available, determine the number threads per process automatically. #if _OPENMP SET_INT_PROP(FlowModelParameters, ThreadsPerProcess, -1); @@ -111,6 +120,12 @@ namespace Opm double tolerance_well_control_; /// Tolerance for the pressure equations for multisegment wells double tolerance_pressure_ms_wells_; + + /// Relaxed tolerance for the inner iteration for the MSW flow solution + double relaxed_inner_tolerance_flow_ms_well_; + /// Relaxed tolerance for the inner iteration for the MSW pressure solution + double relaxed_inner_tolerance_pressure_ms_well_; + /// Maximum pressure change over an iteratio for ms wells double max_pressure_change_ms_wells_; @@ -120,6 +135,12 @@ namespace Opm /// Maximum inner iteration number for ms wells int max_inner_iter_ms_wells_; + /// Strict inner iteration number for ms wells + int strict_inner_iter_ms_wells_; + + /// Regularization factor for ms wells + int regularization_factor_ms_wells_; + /// Maximum iteration number of the well equation solution int max_welleq_iter_; @@ -166,9 +187,13 @@ namespace Opm max_welleq_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxWelleqIter); use_multisegment_well_ = EWOMS_GET_PARAM(TypeTag, bool, UseMultisegmentWell); tolerance_pressure_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, TolerancePressureMsWells); + relaxed_inner_tolerance_flow_ms_well_ = EWOMS_GET_PARAM(TypeTag, Scalar, RelaxedFlowTolInnerIterMsw); + relaxed_inner_tolerance_pressure_ms_well_ = EWOMS_GET_PARAM(TypeTag, Scalar, RelaxedPressureTolInnerIterMsw); max_pressure_change_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells); use_inner_iterations_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseInnerIterationsMsWells); max_inner_iter_ms_wells_ = EWOMS_GET_PARAM(TypeTag, int, MaxInnerIterMsWells); + strict_inner_iter_ms_wells_ = EWOMS_GET_PARAM(TypeTag, int, StrictInnerIterMsWells); + regularization_factor_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, RegularizationFactorMsw); maxSinglePrecisionTimeStep_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays) *24*60*60; max_strict_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxStrictIter); solve_welleq_initially_ = EWOMS_GET_PARAM(TypeTag, bool, SolveWelleqInitially); @@ -192,9 +217,13 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, int, MaxWelleqIter, "Maximum number of iterations to determine solution the well equations"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseMultisegmentWell, "Use the well model for multi-segment wells instead of the one for single-segment wells"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, TolerancePressureMsWells, "Tolerance for the pressure equations for multi-segment wells"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, RelaxedFlowTolInnerIterMsw, "Relaxed tolerance for the inner iteration for the MSW flow solution"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, RelaxedPressureTolInnerIterMsw, "Relaxed tolerance for the inner iteration for the MSW pressure solution"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells, "Maximum relative pressure change for a single iteration of the multi-segment well model"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseInnerIterationsMsWells, "Use nested iterations for multi-segment wells"); EWOMS_REGISTER_PARAM(TypeTag, int, MaxInnerIterMsWells, "Maximum number of inner iterations for multi-segment wells"); + EWOMS_REGISTER_PARAM(TypeTag, int, StrictInnerIterMsWells, "Number of inner iterations for multi-segment wells with strict tolerance"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, RegularizationFactorMsw, "Regularization factor for ms wells"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays, "Maximum time step size where single precision floating point arithmetic can be used solving for the linear systems of equations"); EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIter, "Maximum number of Newton iterations before relaxed tolerances are used for the CNV convergence criterion"); EWOMS_REGISTER_PARAM(TypeTag, bool, SolveWelleqInitially, "Fully solve the well equations before each iteration of the reservoir model"); diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 9cc9b12f9..f24bcd14b 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -128,7 +128,7 @@ namespace Opm Opm::DeferredLogger& deferred_logger) const override; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) const override; + virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance = false) const override; /// Ax = Ax - C D^-1 B x virtual void apply(const BVector& x, BVector& Ax) const override; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index e85ba1085..21f547222 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -536,7 +536,7 @@ namespace Opm template ConvergenceReport MultisegmentWell:: - getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) const + getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance) const { assert(int(B_avg.size()) == num_components_); @@ -577,11 +577,14 @@ namespace Opm if (eq_idx < num_components_) { // phase or component mass equations const double flux_residual = maximum_residual[eq_idx]; // TODO: the report can not handle the segment number yet. + if (std::isnan(flux_residual)) { report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, name()}); } else if (flux_residual > param_.max_residual_allowed_) { report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, name()}); - } else if (flux_residual > param_.tolerance_wells_) { + } else if (!relax_tolerance && flux_residual > param_.tolerance_wells_) { + report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()}); + } else if (flux_residual > param_.relaxed_inner_tolerance_flow_ms_well_) { report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()}); } } else { // pressure equation @@ -591,7 +594,9 @@ namespace Opm report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()}); } else if (std::isinf(pressure_residual)) { report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, name()}); - } else if (pressure_residual > param_.tolerance_pressure_ms_wells_) { + } else if (!relax_tolerance && pressure_residual > param_.tolerance_pressure_ms_wells_) { + report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()}); + } else if (pressure_residual > param_.relaxed_inner_tolerance_pressure_ms_well_) { report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()}); } } @@ -1119,8 +1124,8 @@ namespace Opm // update the segment pressure { const int sign = dwells[seg][SPres] > 0.? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); - // const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); + //const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); + const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); } @@ -2369,17 +2374,20 @@ namespace Opm int it = 0; // relaxation factor double relaxation_factor = 1.; - const double min_relaxation_factor = 0.2; + const double min_relaxation_factor = 0.6; bool converged = false; int stagnate_count = 0; + bool relax_convergence = false; for (; it < max_iter_number; ++it, ++debug_cost_counter_) { assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger); const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + if (it > param_.strict_inner_iter_ms_wells_) + relax_convergence = true; - const auto report = getWellConvergence(well_state, B_avg, deferred_logger); + const auto report = getWellConvergence(well_state, B_avg, deferred_logger, relax_convergence); if (report.converged()) { converged = true; break; @@ -2395,19 +2403,22 @@ namespace Opm // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor, // for example, to recover it to be bigger - if (!is_stagnate) { - stagnate_count = 0; - } if (is_oscillate || is_stagnate) { // HACK! - if (is_stagnate && relaxation_factor == min_relaxation_factor) { + std::ostringstream sstr; + if (relaxation_factor == min_relaxation_factor) { // Still stagnating, terminate iterations if 5 iterations pass. ++stagnate_count; - if (stagnate_count == 5) { - // break; + if (stagnate_count == 6) { + sstr << " well " << name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n"; + const auto reportStag = getWellConvergence(well_state, B_avg, deferred_logger, true); + if (reportStag.converged()) { + converged = true; + sstr << " well " << name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations"; + deferred_logger.debug(sstr.str()); + return; + } } - } else { - stagnate_count = 0; } // a factor value to reduce the relaxation_factor @@ -2415,7 +2426,6 @@ namespace Opm relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor); // debug output - std::ostringstream sstr; if (is_stagnate) { sstr << " well " << name() << " observes stagnation in inner iteration " << it << "\n"; @@ -2433,7 +2443,9 @@ namespace Opm // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state if (converged) { std::ostringstream sstr; - sstr << " well " << name() << " manage to get converged within " << it << " inner iterations"; + sstr << " well " << name() << " manage to get converged within " << it << " inner iterations \n"; + if (relax_convergence) + sstr << "A relaxed tolerance is used after "<< param_.strict_inner_iter_ms_wells_ << " iterations"; deferred_logger.debug(sstr.str()); } else { std::ostringstream sstr; @@ -2495,9 +2507,14 @@ namespace Opm // TODO: without considering the efficiencty factor for now { const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg); + + // Add a regularization_factor to increase the accumulation term + // This will make the system less stiff and help convergence for + // difficult cases + const Scalar regularization_factor = param_.regularization_factor_ms_wells_; // for each component for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell accumulation_term = (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx) + const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx) - segment_fluid_initial_[seg][comp_idx]) / dt; resWell_[seg][comp_idx] += accumulation_term.value(); @@ -2818,7 +2835,7 @@ namespace Opm vol_ratio += mix[comp_idx] / b[comp_idx]; } - // segment volume + // We increase the segment volume with a factor 10 to stabilize the system. const double volume = segmentSet()[seg_idx].volume(); return volume / vol_ratio; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index c29ec0593..7ed6c19f4 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -177,7 +177,8 @@ namespace Opm /// check whether the well equations get converged for this well virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, - Opm::DeferredLogger& deferred_logger) const override; + Opm::DeferredLogger& deferred_logger, + const bool relax_tolerance = false) const override; /// Ax = Ax - C D^-1 B x virtual void apply(const BVector& x, BVector& Ax) const override; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index d1c463599..c1a96944b 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2087,7 +2087,8 @@ namespace Opm StandardWell:: getWellConvergence(const WellState& well_state, const std::vector& B_avg, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger, + const bool /*relax_tolerance*/) const { // the following implementation assume that the polymer is always after the w-o-g phases // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index d22ff9d4f..abc915096 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -157,7 +157,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const = 0; - virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger) const = 0; + virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance = false) const = 0; virtual void solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0;