diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 9cd736924..737793b92 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -105,7 +105,7 @@ try Opm::LinearSolverFactory linsolver(param); Opm::NewtonIterationBlackoilSimple fis_solver(linsolver); - Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, fis_solver); + Opm::FullyImplicitBlackoilSolver solver(param, *g, props, geo, 0, *wells, fis_solver); Opm::BlackoilState state; initStateBasic(*g, props0, param, 0.0, state); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 29aebeac7..481df035c 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -31,6 +31,7 @@ struct Wells; namespace Opm { + namespace parameter { class ParameterGroup; } class DerivedGeology; class RockCompressibility; class NewtonIterationBlackoilInterface; @@ -56,13 +57,15 @@ namespace Opm { /// Construct a solver. It will retain references to the /// arguments of this functions, and they are expected to /// remain in scope for the lifetime of the solver. + /// \param[in] param parameters /// \param[in] grid grid data structure /// \param[in] fluid fluid properties /// \param[in] geo rock properties /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] wells well structure /// \param[in] linsolver linear solver - FullyImplicitBlackoilSolver(const Grid& grid , + FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param, + const Grid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, @@ -141,6 +144,9 @@ namespace Opm { HelperOps ops_; const WellOps wops_; const M grav_; + double dp_max_rel_; + double ds_max_; + double drs_max_rel_; std::vector rq_; std::vector phaseCondition_; @@ -278,6 +284,9 @@ namespace Opm { void stablizeNewton(V &dx, V &dxOld, const bool &oscillate, const bool &stagnate, const double omega, const RelaxType relax_type) const; + const double dpMaxRel () const { return dp_max_rel_; } + const double dsMax () const { return ds_max_; } + const double drsMaxRel () const { return drs_max_rel_; } }; } // namespace Opm diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 805ac8e1b..2d950b060 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -212,7 +213,8 @@ namespace { template FullyImplicitBlackoilSolver:: - FullyImplicitBlackoilSolver(const Grid& grid , + FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param, + const Grid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, @@ -230,17 +232,27 @@ namespace { , ops_ (grid) , wops_ (wells) , grav_ (gravityOperator(grid_, ops_, geo_)) + , dp_max_rel_ ( 1.0e9 ) + , ds_max_ ( 0.2 ) + , drs_max_rel_ ( 1.0e9 ) , rq_ (fluid.numPhases()) , phaseCondition_(AutoDiffGrid::numCells(grid)) , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), ADB::null() } ) { + if (param.has("dp_max_rel")) { + dp_max_rel_ = param.get(std::string("dp_max_rel")); + } + if (param.has("ds_max")) { + ds_max_ = param.get("ds_max"); + } + if (param.has("drs_max_rel")) { + drs_max_rel_ = param.get("drs_max_rel"); + } } - - template void FullyImplicitBlackoilSolver:: @@ -1251,7 +1263,7 @@ namespace { assert(varstart == dx.size()); // Pressure update. - const double dpmaxrel = 0.8; + const double dpmaxrel = dpMaxRel(); const V p_old = Eigen::Map(&state.pressure()[0], nc, 1); const V absdpmax = dpmaxrel*p_old.abs(); const V dp_limited = sign(dp) * dp.abs().min(absdpmax); @@ -1262,40 +1274,53 @@ namespace { // Saturation updates. const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np); - const double dsmax = 0.3; + const double dsmax = dsMax(); V so = one; V sw; + V sg; - if (active_[ Water ]) { - const int pos = pu.phase_pos[ Water ]; - const V sw_old = s_old.col(pos); - const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax); - sw = (sw_old - dsw_limited).unaryExpr(Chop01()); - so -= sw; - for (int c = 0; c < nc; ++c) { - state.saturation()[c*np + pos] = sw[c]; + { + V maxVal = zero; + V dso = zero; + if (active_[Water]){ + maxVal = dsw.abs().max(maxVal); + dso = dso - dsw; + } + + V dsg; + if (active_[Gas]){ + dsg = isSg * dxvar - isRv * dsw; + maxVal = dsg.abs().max(maxVal); + dso = dso - dsg; + } + + maxVal = dso.abs().max(maxVal); + + V step = dsmax/maxVal; + step = step.min(1.); + + if (active_[Water]) { + const int pos = pu.phase_pos[ Water ]; + const V sw_old = s_old.col(pos); + sw = sw_old - step * dsw; + so -= sw; + } + + if (active_[Gas]) { + const int pos = pu.phase_pos[ Gas ]; + const V sg_old = s_old.col(pos); + sg = sg_old - step * dsg; + so -= sg; } } - V sg; - if (active_[Gas]) { - const int pos = pu.phase_pos[ Gas ]; - const V sg_old = s_old.col(pos); - const V dsg = isSg * dxvar - isRv * dsw; - const V dsg_limited = sign(dsg) * dsg.abs().min(dsmax); - sg = sg_old - dsg_limited; - so -= sg; - } - - - - const double drsmax = 1e9; + const double drsmaxrel = drsMaxRel(); const double drvmax = 1e9;//% same as in Mrst V rs; if (disgas) { const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); const V drs = isRs * dxvar; - const V drs_limited = sign(drs) * drs.abs().min(drsmax); + const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel); rs = rs_old - drs_limited; } V rv; @@ -1431,7 +1456,8 @@ namespace { // Bhp update. const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V bhp = bhp_old - dbhp; + const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + const V bhp = bhp_old - dbhp_limited; std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index d40d2e0cd..81fa58fce 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -200,7 +200,7 @@ namespace Opm wells_(wells_manager.c_wells()), gravity_(gravity), geo_(grid_, props_, gravity_), - solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) + solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) /* param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10),