From e9e4416428fefec92242429978995f8c78f5bdb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 4 Aug 2017 11:39:09 +0200 Subject: [PATCH 1/7] Make computeCellState() a const method. --- opm/autodiff/BlackoilReorderingTransportModel.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index 700822536..3764c9df1 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -460,7 +460,7 @@ namespace Opm { template - void computeCellState(const int cell, const State& state, CellState& cstate) + void computeCellState(const int cell, const State& state, CellState& cstate) const { assert(numPhases() == 3); // I apologize for this to my future self, that will have to fix it. From c39341006de9b708e967ee0438ec6a937f967db9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 4 Aug 2017 11:39:51 +0200 Subject: [PATCH 2/7] Bugfix: follow opm-material convention for capillary pressure --- opm/autodiff/BlackoilReorderingTransportModel.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index 3764c9df1..0f4645ebb 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -499,8 +499,8 @@ namespace Opm { // Compute phase pressures. cstate.p[Oil] = constant(poval); - cstate.p[Water] = cstate.p[Oil] - cstate.pc[Water]; // pcow = po - pw - cstate.p[Gas] = cstate.p[Oil] + cstate.pc[Gas]; // pcog = pg - po (!) + cstate.p[Water] = cstate.p[Oil] + cstate.pc[Water]; // pcow = pw - po (!) [different from old convention] + cstate.p[Gas] = cstate.p[Oil] + cstate.pc[Gas]; // pcog = pg - po // Compute PVT properties. cstate.temperature = constant(0.0); // Temperature is not used. From 608a674858693c27e56ee3cab3b60c2092262c77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 7 Aug 2017 16:17:57 +0200 Subject: [PATCH 3/7] Make dh_sat behave like in coupled transport solver. --- opm/autodiff/BlackoilReorderingTransportModel.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index 0f4645ebb..0ba036a9f 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -753,14 +753,15 @@ namespace Opm { // values from cstate_[other]. Eval dh[3]; Eval dh_sat[3]; + const Eval grad_oil_press = cstate_[other].p[Oil] - st.p[Oil]; for (int phase : { Water, Oil, Gas }) { const Eval gradp = cstate_[other].p[phase] - st.p[phase]; const Eval rhoavg = 0.5 * (st.rho[phase] + cstate_[other].rho[phase]); dh[phase] = gradp - rhoavg * gdz; - dh_sat[phase] = rhoavg * gdz; // TODO: not equivalent to formulation in BlackoilTransportModel for cap. press. if (Base::use_threshold_pressure_) { - applyThresholdPressure(conn.index, dh[phase]); // TODO: Should also dh_sat be treated here? Also: dh is unused, this has no effect! + applyThresholdPressure(conn.index, dh[phase]); } + dh_sat[phase] = grad_oil_press - dh[phase]; } const double tran = trans_all_[conn.index]; // TODO: include tr_mult effect. const auto& m1 = st.lambda; From 1e0facec4ad0147c59ca579be514170904f98946 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Aug 2017 15:47:39 +0200 Subject: [PATCH 4/7] Fix updateState() saturation update. --- opm/autodiff/BlackoilReorderingTransportModel.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index 0ba036a9f..2884b98d9 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -855,13 +855,14 @@ namespace Opm { // Get saturation updates. const double dsw = dx[0]; - double dso = -dsw; double dsg = 0.0; auto& hcstate = state_.reservoir_state.hydroCarbonState()[cell]; if (hcstate == HydroCarbonState::GasAndOil) { dsg = dx[1]; - dso -= dsg; + } else if (hcstate == HydroCarbonState::GasOnly) { + dsg = -dsw; } + const double dso = -(dsw + dsg); // Handle too large saturation changes. const double maxval = std::max(std::fabs(dsw), std::max(std::fabs(dso), std::fabs(dsg))); From 9cf68321405fcd387c006bca1d987f146e8a371e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Aug 2017 11:25:47 +0200 Subject: [PATCH 5/7] Ignore limits to rs and rv changes. As implemented with a relative limit, even with 1e9 default limit it would still be impossible to get away from a zero value. It is possible that the limits may return later, implemented in a less buggy way, however for now they do not seem necessary. --- .../BlackoilReorderingTransportModel.hpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index 2884b98d9..a34ad0630 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -888,20 +888,24 @@ namespace Opm { double& rs = state_.reservoir_state.gasoilratio()[cell]; const double rs_old = rs; if (hcstate == HydroCarbonState::OilOnly) { - const double max_allowed_change = std::fabs(rs_old) * Base::drMaxRel(); + // const double max_allowed_change = std::fabs(rs_old) * Base::drMaxRel(); const double drs = dx[1]; - const double factor = std::min(1.0, max_allowed_change / std::fabs(drs)); - rs += factor*drs; + // const double factor = std::min(1.0, max_allowed_change / std::fabs(drs)); + // rs += factor*drs; + rs += drs; + rs = std::max(rs, 0.0); } // Update rv. double& rv = state_.reservoir_state.rv()[cell]; const double rv_old = rv; if (hcstate == HydroCarbonState::GasOnly) { - const double max_allowed_change = std::fabs(rv_old) * Base::drMaxRel(); + // const double max_allowed_change = std::fabs(rv_old) * Base::drMaxRel(); const double drv = dx[1]; - const double factor = std::min(1.0, max_allowed_change / std::fabs(drv)); - rv += factor*drv; + // const double factor = std::min(1.0, max_allowed_change / std::fabs(drv)); + // rv += factor*drv; + rv += drv; + rv = std::max(rv, 0.0); } const double epsilon = std::sqrt(std::numeric_limits::epsilon()); From 784bcf6892d8880eb97a084e2f9db98e51fa232f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Aug 2017 11:40:05 +0200 Subject: [PATCH 6/7] Experiment with relaxing single cell iterations. --- .../BlackoilReorderingTransportModel.hpp | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index a34ad0630..f9b1efee2 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -642,23 +642,37 @@ namespace Opm { // Newton loop. int iter = 0; - const int max_iter = 25; + const int max_iter = 100; double relaxation = 1.0; while (!getConvergence(cell, res) && iter < max_iter) { Vec2 dx; jac.solve(dx, res); dx *= relaxation; + const auto hcstate_old = state_.reservoir_state.hydroCarbonState()[cell]; updateState(cell, -dx); + const auto hcstate = state_.reservoir_state.hydroCarbonState()[cell]; assembleSingleCell(cell, res, jac); ++iter; - // if (iter > 15) { - // relaxation = 0.7; - // std::ostringstream os; - // os << "Iteration " << iter << " in cell " << cell << ", residual = " << res - // << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas] - // << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << "}, dx = " << dx; - // OpmLog::debug(os.str()); - // } + if (iter > 10) { + relaxation = 0.85; + if (iter > 15) { + relaxation = 0.70; + } + if (iter > 20) { + relaxation = 0.55; + } + if (iter > 25) { + relaxation = 0.40; + } + if (iter > 30) { + relaxation = 0.25; + } + std::ostringstream os; + os << "Iteration " << iter << " in cell " << cell << ", residual = " << res + << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas] + << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << "}, dx = " << dx << ", hcstate: " << hcstate_old << " -> " << hcstate; + OpmLog::debug(os.str()); + } } if (iter == max_iter) { std::ostringstream os; From 67823c5893301d1b83a4884b572482116654d5ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Aug 2017 13:09:22 +0200 Subject: [PATCH 7/7] Disable debugging output of detailed cell data. --- opm/autodiff/BlackoilReorderingTransportModel.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index f9b1efee2..cfd320c7d 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -667,11 +667,11 @@ namespace Opm { if (iter > 30) { relaxation = 0.25; } - std::ostringstream os; - os << "Iteration " << iter << " in cell " << cell << ", residual = " << res - << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas] - << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << "}, dx = " << dx << ", hcstate: " << hcstate_old << " -> " << hcstate; - OpmLog::debug(os.str()); + // std::ostringstream os; + // os << "Iteration " << iter << " in cell " << cell << ", residual = " << res + // << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas] + // << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << "}, dx = " << dx << ", hcstate: " << hcstate_old << " -> " << hcstate; + // OpmLog::debug(os.str()); } } if (iter == max_iter) {