From 4cef1510915a51431f546933c485f3560715b4ba Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 13:26:44 +0200 Subject: [PATCH 1/9] Add dsmax for the update of saturations. dsmax is the absolute limit for saturation update. --- .../FullyImplicitBlackoilSolver_impl.hpp | 55 ++++++++++++------- 1 file changed, 34 insertions(+), 21 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 24e57e925..5d3b8e212 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1235,33 +1235,46 @@ 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 = 0.05; 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 drvmax = 1e9;//% same as in Mrst V rs; From b001c580b94104d80d879a1443a6fcb778b939aa Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 13:32:34 +0200 Subject: [PATCH 2/9] Add drsmaxrel for the update of rs. drsmaxrel is a relative limit. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 5d3b8e212..3d15f5994 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1275,13 +1275,13 @@ namespace { } } - const double drsmax = 1e9; + const double drsmaxrel = 0.2; 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; From 55b0164d4ad895617902e7cf392b4600804bdacc Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 13:46:00 +0200 Subject: [PATCH 3/9] Apply dpmaxrel to the update of bhp And also change the default value for dpmaxrel. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 3d15f5994..6e330ccad 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1224,7 +1224,7 @@ namespace { assert(varstart == dx.size()); // Pressure update. - const double dpmaxrel = 0.8; + const double dpmaxrel = 0.2; 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); @@ -1417,7 +1417,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()); } From 94e3fd3fcb416fbd54c468ee61a1048eb1f1cb89 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 13:52:33 +0200 Subject: [PATCH 4/9] Add more residual output in getConvergence. Output the residualWellFlux, residualWell, MB to monitor the convergence process. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 6e330ccad..ebdf6e0d2 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1682,11 +1682,13 @@ namespace { bool converged = converged_MB && converged_CNV && converged_Well; -#ifdef OPM_VERBOSE +// #ifdef OPM_VERBOSE + std::cout << " residualWellFlux " << residualWellFlux << " residualWell " << residualWell << std::endl; std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl; + std::cout << " MB " << fabs(BW_avg*RW_sum) << " " << fabs(BO_avg*RO_sum) << " " << fabs(BG_avg*RG_sum) << std::endl; std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV << " converged_Well " << converged_Well << " converged " << converged << std::endl; -#endif +// #endif return converged; } From cd50b54ddf59bc4e9f4ceab7fa001c36afaa7b02 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 18:02:55 +0200 Subject: [PATCH 5/9] Finishing the modification for solver class. --- examples/test_implicit_ad.cpp | 2 +- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 12 ++++++++- .../FullyImplicitBlackoilSolver_impl.hpp | 25 +++++++++++++------ .../SimulatorFullyImplicitBlackoil_impl.hpp | 2 +- 4 files changed, 31 insertions(+), 10 deletions(-) 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 d3afee547..e4f26af08 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,16 @@ 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, @@ -137,6 +141,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_; @@ -266,6 +273,9 @@ namespace Opm { /// residual mass balance (tol_cnv). bool getConvergence(const double dt); + 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 ebdf6e0d2..52983a7d7 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -209,10 +210,10 @@ namespace { } // Anonymous 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, @@ -235,12 +236,22 @@ namespace { , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), ADB::null() } ) + , dp_max_rel_ ( 0.2 ) + , ds_max_ ( 0.05 ) + , drs_max_rel_ ( 0.2 ) { + if (param.has("dp_max_rel")) { + dp_max_rel_ = param.get(std::string("dp_max_rel")); + } + if (param.has("dp_max")) { + ds_max_ = param.get("dp_max"); + } + if (param.has("drs_max_rel")) { + ds_max_ = param.get("drs_max_rel"); + } } - - template void FullyImplicitBlackoilSolver:: @@ -1224,7 +1235,7 @@ namespace { assert(varstart == dx.size()); // Pressure update. - const double dpmaxrel = 0.2; + 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); @@ -1235,7 +1246,7 @@ namespace { // Saturation updates. const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np); - const double dsmax = 0.05; + const double dsmax = dsMax(); V so = one; V sw; V sg; @@ -1275,7 +1286,7 @@ namespace { } } - const double drsmaxrel = 0.2; + const double drsmaxrel = drsMaxRel(); const double drvmax = 1e9;//% same as in Mrst V rs; if (disgas) { 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), From f2ecbf163e216c7f3bf25a1621eef8042106089d Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 18:13:36 +0200 Subject: [PATCH 6/9] Cleaning up some debugging output. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 52983a7d7..e649f2f47 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -243,11 +243,11 @@ namespace { if (param.has("dp_max_rel")) { dp_max_rel_ = param.get(std::string("dp_max_rel")); } - if (param.has("dp_max")) { - ds_max_ = param.get("dp_max"); + if (param.has("ds_max")) { + ds_max_ = param.get("ds_max"); } if (param.has("drs_max_rel")) { - ds_max_ = param.get("drs_max_rel"); + drs_max_rel_ = param.get("drs_max_rel"); } } @@ -1693,13 +1693,13 @@ namespace { bool converged = converged_MB && converged_CNV && converged_Well; -// #ifdef OPM_VERBOSE +#ifdef OPM_VERBOSE std::cout << " residualWellFlux " << residualWellFlux << " residualWell " << residualWell << std::endl; std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl; std::cout << " MB " << fabs(BW_avg*RW_sum) << " " << fabs(BO_avg*RO_sum) << " " << fabs(BG_avg*RG_sum) << std::endl; std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV << " converged_Well " << converged_Well << " converged " << converged << std::endl; -// #endif +#endif return converged; } From d461eb9e761fa68691de0c4650ecb78573aad4f1 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 18:15:42 +0200 Subject: [PATCH 7/9] Cleaning up some debugging output. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index e649f2f47..cf392a4ee 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1694,9 +1694,7 @@ namespace { bool converged = converged_MB && converged_CNV && converged_Well; #ifdef OPM_VERBOSE - std::cout << " residualWellFlux " << residualWellFlux << " residualWell " << residualWell << std::endl; std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl; - std::cout << " MB " << fabs(BW_avg*RW_sum) << " " << fabs(BO_avg*RO_sum) << " " << fabs(BG_avg*RG_sum) << std::endl; std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV << " converged_Well " << converged_Well << " converged " << converged << std::endl; #endif From 2b31ab6111e70460b365dd3ab405e6fb63359eec Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 16 May 2014 18:27:23 +0200 Subject: [PATCH 8/9] Removing one blank line. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 1 - opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index e4f26af08..df1126eed 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -64,7 +64,6 @@ 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 parameter::ParameterGroup& param, const Grid& grid , const BlackoilPropsAdInterface& fluid, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index cf392a4ee..195ff4d9f 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -210,6 +210,7 @@ namespace { } // Anonymous namespace + template FullyImplicitBlackoilSolver:: FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param, From 79916078f285337fad3176eb0b7a5f9df232ea77 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 20 May 2014 13:05:11 +0200 Subject: [PATCH 9/9] Changing the defalut value for the paramters. Add reordering the intialization order to removing the reordered warning. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 195ff4d9f..e2bfa67f3 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -232,14 +232,14 @@ 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() } ) - , dp_max_rel_ ( 0.2 ) - , ds_max_ ( 0.05 ) - , drs_max_rel_ ( 0.2 ) { if (param.has("dp_max_rel")) { dp_max_rel_ = param.get(std::string("dp_max_rel"));