From 13d3a4244814dd0862fc81ac276bdbe931df4b04 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 13:39:11 +0200 Subject: [PATCH 1/6] Fix phase property references Methods 'computeRelperm()' and 'computePressures()' *always* return a three-element vector of phase properties. We must therefore translate to canonical phase indices before indexing into the results. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 4ac2cb7e2..814145c7c 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -711,7 +711,7 @@ namespace { const std::vector kr = computeRelPerm(state); const std::vector pressures = computePressures(state); for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { - computeMassFlux(phaseIdx, transi, kr[phaseIdx], pressures[phaseIdx], state); + computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], pressures[canph_[phaseIdx]], state); // std::cout << "===== kr[" << phase << "] = \n" << std::endl; // std::cout << kr[phase]; // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl; From 35d883319b25100270f5889c8bfabf0751b42ff4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 16:49:57 +0200 Subject: [PATCH 2/6] variableState: Don't reference non-existing variable The "xvar" exists only if Gas is active. Therefore, we cannot extract that variable from "vars" unless we know that Gas is an active phase. Failing to do so would wrongfully increment 'nextvar' whence the final BHP variable would be an out-of-bounds access. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 814145c7c..670b5f43a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -532,8 +532,8 @@ namespace { std::vector all_cells = buildAllCells(nc); ADB rsSat = fluidRsSat(state.pressure,all_cells); ADB rvSat = fluidRvSat(state.pressure,all_cells); - ADB xvar = vars[ nextvar++ ]; if (active_[ Gas]) { + ADB xvar = vars[ nextvar++ ]; ADB sg = isSg*xvar + isRv* so; state.saturation[ pu.phase_pos[ Gas ] ] = sg; so = so - sg; From f69530a7d41d2ab8473caea0875c865e96825eaf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 16:51:54 +0200 Subject: [PATCH 3/6] variableState: Defer r{s,v}Sat calculation Quantities rsSat and rvSat are not needed unless we have an active Gas phase. Don't calculate them until we know that they are actually needed. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 670b5f43a..d439ad2c5 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -528,11 +528,11 @@ namespace { so = so - sw; } - // Define Sg Rs and Rv in terms of xvar. - std::vector all_cells = buildAllCells(nc); - ADB rsSat = fluidRsSat(state.pressure,all_cells); - ADB rvSat = fluidRvSat(state.pressure,all_cells); if (active_[ Gas]) { + // Define Sg Rs and Rv in terms of xvar. + std::vector all_cells = buildAllCells(nc); + ADB rsSat = fluidRsSat(state.pressure,all_cells); + ADB rvSat = fluidRvSat(state.pressure,all_cells); ADB xvar = vars[ nextvar++ ]; ADB sg = isSg*xvar + isRv* so; state.saturation[ pu.phase_pos[ Gas ] ] = sg; From eccc5c96946af48daeb568507534b9b3740bcd95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 17:22:50 +0200 Subject: [PATCH 4/6] variableState: Use pre-computed list of "all cells" Class FullyImplicitBlackoilSolver already features a list of "all" cells, built at object construction time. There's no need to re-compute that list on every call to variableState() (when Gas is active). --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index d439ad2c5..12cb372ef 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -530,9 +530,8 @@ namespace { if (active_[ Gas]) { // Define Sg Rs and Rv in terms of xvar. - std::vector all_cells = buildAllCells(nc); - ADB rsSat = fluidRsSat(state.pressure,all_cells); - ADB rvSat = fluidRvSat(state.pressure,all_cells); + ADB rsSat = fluidRsSat(state.pressure, cells_); + ADB rvSat = fluidRvSat(state.pressure, cells_); ADB xvar = vars[ nextvar++ ]; ADB sg = isSg*xvar + isRv* so; state.saturation[ pu.phase_pos[ Gas ] ] = sg; From 78b87b88573b07189d328617c7c870dbfcff2b6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 18:42:58 +0200 Subject: [PATCH 5/6] variableState: Don't form objects when references suffice There's no need to form new objects, especially for the 'xvar', when there are existing objects to which read-only references can be bound. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 12cb372ef..52bd1ec2c 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -530,10 +530,10 @@ namespace { if (active_[ Gas]) { // Define Sg Rs and Rv in terms of xvar. - ADB rsSat = fluidRsSat(state.pressure, cells_); - ADB rvSat = fluidRvSat(state.pressure, cells_); - ADB xvar = vars[ nextvar++ ]; - ADB sg = isSg*xvar + isRv* so; + const ADB& rsSat = fluidRsSat(state.pressure, cells_); + const ADB& rvSat = fluidRvSat(state.pressure, cells_); + const ADB& xvar = vars[ nextvar++ ]; + const ADB& sg = isSg*xvar + isRv* so; state.saturation[ pu.phase_pos[ Gas ] ] = sg; so = so - sg; From a5d30170628d30dd04392e9be44317b3855a034f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 27 May 2014 18:54:17 +0200 Subject: [PATCH 6/6] assemble: Fix out-of-bounds indexing The 'material_balance_eq' is indexed by active, not canonical, phase indices. Replace 'Oil' and 'Gas' with the appropriate 'phase_pos' values. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 52bd1ec2c..88c03bb7a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -732,18 +732,18 @@ namespace { // The extra terms in the accumulation part of the equation are already handled. if (active_[ Oil ] && active_[ Gas ]) { const int po = fluid_.phaseUsage().phase_pos[ Oil ]; + const int pg = fluid_.phaseUsage().phase_pos[ Gas ]; + const UpwindSelector upwindOil(grid_, ops_, rq_[po].head.value()); const ADB rs_face = upwindOil.select(state.rs); - residual_.material_balance_eq[ Gas ] += ops_.div * (rs_face * rq_[po].mflux); - - const int pg = fluid_.phaseUsage().phase_pos[ Gas ]; const UpwindSelector upwindGas(grid_, ops_, rq_[pg].head.value()); const ADB rv_face = upwindGas.select(state.rv); - residual_.material_balance_eq[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux); + residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux); + residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux); // DUMP(residual_.material_balance_eq[ Gas ]);