From 80c38d5a1a9a95179d62a997d8fbe012d1caa19d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 7 Jul 2016 13:04:57 +0200 Subject: [PATCH 1/5] Add connectionMultiPhaseUpwind() free function. --- CMakeLists_files.cmake | 2 + opm/autodiff/multiPhaseUpwind.cpp | 87 +++++++++++++++++++++++++++++++ opm/autodiff/multiPhaseUpwind.hpp | 45 ++++++++++++++++ 3 files changed, 134 insertions(+) create mode 100644 opm/autodiff/multiPhaseUpwind.cpp create mode 100644 opm/autodiff/multiPhaseUpwind.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2392fa88d..bf91c3d9b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -35,6 +35,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/GridHelpers.cpp opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/moduleVersion.cpp + opm/autodiff/multiPhaseUpwind.cpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp opm/autodiff/SimulatorIncompTwophaseAd.cpp opm/autodiff/TransportSolverTwophaseAd.cpp @@ -188,6 +189,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/ISTLSolver.hpp opm/autodiff/IterationReport.hpp opm/autodiff/moduleVersion.hpp + opm/autodiff/multiPhaseUpwind.hpp opm/autodiff/NewtonIterationBlackoilCPR.hpp opm/autodiff/NewtonIterationBlackoilInterface.hpp opm/autodiff/NewtonIterationBlackoilInterleaved.hpp diff --git a/opm/autodiff/multiPhaseUpwind.cpp b/opm/autodiff/multiPhaseUpwind.cpp new file mode 100644 index 000000000..934208f50 --- /dev/null +++ b/opm/autodiff/multiPhaseUpwind.cpp @@ -0,0 +1,87 @@ +/* + Copyright 2015, 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#include +#include + + +namespace Opm +{ + + + std::array connectionMultiPhaseUpwind(const std::array& head_diff, + const std::array& mob1, + const std::array& mob2, + const double transmissibility, + const double flux) + { + // Based on the paper "Upstream Differencing for Multiphase Flow in Reservoir Simulation", + // by Yann Brenier and Jérôme Jaffré, + // SIAM J. Numer. Anal., 28(3), 685–696. + // DOI:10.1137/0728036 + // + // Notation is based on this paper, except q -> flux, t -> transmissibility. + + enum { NumPhases = 3 }; // TODO: remove this restriction. + + // Get and sort the g values (also called "weights" in the paper) for this connection. + using ValueAndIndex = std::pair; + std::array g; + for (int phase_idx = 0; phase_idx < NumPhases; ++phase_idx) { + g[phase_idx] = ValueAndIndex(head_diff[phase_idx], phase_idx); + } + std::sort(g.begin(), g.end()); + + // Compute theta and r. + // Paper notation: subscript l -> ell (for read/searchability) + // Note that since we index phases from 0, r is one less than in the paper. + std::array theta; + int r = -1; + for (int ell = 0; ell < NumPhases; ++ell) { + theta[ell] = flux; + for (int j = 0; j < NumPhases; ++j) { + if (j < ell) { + theta[ell] += transmissibility * (g[ell].first - g[j].first) * mob2[g[j].second]; + } + if (j > ell) { + theta[ell] += transmissibility * (g[ell].first - g[j].first) * mob1[g[j].second]; + } + } + if (theta[ell] <= 0.0) { + r = ell; + } else { + break; // r is correct, no need to continue + } + } + + // Set upwind array and return. + std::array upwind; + for (int ell = 0; ell < NumPhases; ++ell) { + const int phase_idx = g[ell].second; + upwind[phase_idx] = ell > r ? 1.0 : -1.0; + } + return upwind; + + } + + +} // namespace Opm + diff --git a/opm/autodiff/multiPhaseUpwind.hpp b/opm/autodiff/multiPhaseUpwind.hpp new file mode 100644 index 000000000..8f52824bd --- /dev/null +++ b/opm/autodiff/multiPhaseUpwind.hpp @@ -0,0 +1,45 @@ +/* + Copyright 2015, 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil AS. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_MULTIPHASEUPWIND_HEADER_INCLUDED +#define OPM_MULTIPHASEUPWIND_HEADER_INCLUDED + +#include + +namespace Opm +{ + /// Compute upwind directions for three-phase flow across a connection. + /// + /// @param[in] head_diff head differences by phase + /// @param[in] mob1 phase mobilities for first cell + /// @param[in] mob2 phase mobilities for second cell + /// @param[in] transmissibility tranmissibility of connection + /// @param[in] flux total volume flux across connection + /// @return array containing, for each phase, 1.0 if flow in the + /// direction of the connection, -1.0 if flow in the opposite + /// direction. + std::array connectionMultiPhaseUpwind(const std::array& head_diff, + const std::array& mob1, + const std::array& mob2, + const double transmissibility, + const double flux); +} // namespace Opm + +#endif // OPM_MULTIPHASEUPWIND_HEADER_INCLUDED From f707abae364f922f2fedf24da9d98283198f09c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 1 Jun 2017 13:55:11 +0200 Subject: [PATCH 2/5] Add unit test for connectionMultiPhaseUpwind(). --- CMakeLists_files.cmake | 1 + tests/test_multiphaseupwind.cpp | 107 ++++++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+) create mode 100644 tests/test_multiphaseupwind.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index bf91c3d9b..1d11e874f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -92,6 +92,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_singlecellsolves.cpp tests/test_solventprops_ad.cpp tests/test_multisegmentwells.cpp + tests/test_multiphaseupwind.cpp # tests/test_thresholdpressure.cpp tests/test_wellswitchlogger.cpp tests/test_timer.cpp diff --git a/tests/test_multiphaseupwind.cpp b/tests/test_multiphaseupwind.cpp new file mode 100644 index 000000000..1d9cf60a6 --- /dev/null +++ b/tests/test_multiphaseupwind.cpp @@ -0,0 +1,107 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE OPM-multiPhaseUpwind +#include + +#include + +// For all the following cases we test a setup with two cells, +// forming a gravity column: +// +// ------------------- +// | | +// | Cell 1 | +// | | +// | | | +// ------------------- | flux +// | | V +// | Cell 2 | +// | | +// | | +// ------------------- +// +// The gravity-related head differences are (4, -1, -2) for (w, o, g). +// The mobilities are all 1 and the transmissibility is 1. +// The total flux from cell 1 to 2 will vary from case to case. + +BOOST_AUTO_TEST_CASE(GravityColumnLowFlux) +{ + // Case 1: a gravity column with two cells and low total flux. + // The total flux from cell 1 to 2 is 1.0. + + const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; + const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; + const double transmissibility = 1.0; + const double flux = 1.0; + + const std::array expected_upw = {{ 1.0, -1.0, -1.0 }}; + const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); + BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); + BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); +} + + +BOOST_AUTO_TEST_CASE(GravityColumnMediumFlux) +{ + // Case 2: a gravity column with two cells and medium-sized total flux. + // The total flux from cell 1 to 2 is 5.0. + + const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; + const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; + const double transmissibility = 1.0; + const double flux = 5.0; + + const std::array expected_upw = {{ 1.0, 1.0, -1.0 }}; + const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); + BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); + BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); +} + + +BOOST_AUTO_TEST_CASE(GravityColumnHighFlux) +{ + // Case 3: a gravity column with two cell and high total flux. + // The total flux from cell 1 to 2 is 10.0. + + const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; + const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; + const double transmissibility = 1.0; + const double flux = 10.0; + + const std::array expected_upw = {{ 1.0, 1.0, 1.0 }}; + const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); + BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); + BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); +} From 33025d752eab6a6144d6f647c63c1ebf1b47937f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 1 Jun 2017 14:07:57 +0200 Subject: [PATCH 3/5] Clarified notation in test. --- tests/test_multiphaseupwind.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tests/test_multiphaseupwind.cpp b/tests/test_multiphaseupwind.cpp index 1d9cf60a6..cac7fc806 100644 --- a/tests/test_multiphaseupwind.cpp +++ b/tests/test_multiphaseupwind.cpp @@ -46,7 +46,8 @@ // | | // ------------------- // -// The gravity-related head differences are (4, -1, -2) for (w, o, g). +// The gravity-related head differences gd ~= rho * g * grad z +// are set to (4, -1, -2) for (w, o, g). // The mobilities are all 1 and the transmissibility is 1. // The total flux from cell 1 to 2 will vary from case to case. @@ -55,14 +56,14 @@ BOOST_AUTO_TEST_CASE(GravityColumnLowFlux) // Case 1: a gravity column with two cells and low total flux. // The total flux from cell 1 to 2 is 1.0. - const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array gd = {{ 4.0, -1.0, -2.0 }}; const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; const double transmissibility = 1.0; const double flux = 1.0; const std::array expected_upw = {{ 1.0, -1.0, -1.0 }}; - const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + const std::array upw = Opm::connectionMultiPhaseUpwind(gd, mob1, mob2, transmissibility, flux); BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); @@ -74,14 +75,14 @@ BOOST_AUTO_TEST_CASE(GravityColumnMediumFlux) // Case 2: a gravity column with two cells and medium-sized total flux. // The total flux from cell 1 to 2 is 5.0. - const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array gd = {{ 4.0, -1.0, -2.0 }}; const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; const double transmissibility = 1.0; const double flux = 5.0; const std::array expected_upw = {{ 1.0, 1.0, -1.0 }}; - const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + const std::array upw = Opm::connectionMultiPhaseUpwind(gd, mob1, mob2, transmissibility, flux); BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); @@ -93,14 +94,14 @@ BOOST_AUTO_TEST_CASE(GravityColumnHighFlux) // Case 3: a gravity column with two cell and high total flux. // The total flux from cell 1 to 2 is 10.0. - const std::array head_diff = {{ 4.0, -1.0, -2.0 }}; + const std::array gd = {{ 4.0, -1.0, -2.0 }}; const std::array mob1 = {{ 1.0, 1.0, 1.0 }}; const std::array mob2 = {{ 1.0, 1.0, 1.0 }}; const double transmissibility = 1.0; const double flux = 10.0; const std::array expected_upw = {{ 1.0, 1.0, 1.0 }}; - const std::array upw = Opm::connectionMultiPhaseUpwind(head_diff, mob1, mob2, transmissibility, flux); + const std::array upw = Opm::connectionMultiPhaseUpwind(gd, mob1, mob2, transmissibility, flux); BOOST_CHECK_EQUAL(upw[0], expected_upw[0]); BOOST_CHECK_EQUAL(upw[1], expected_upw[1]); BOOST_CHECK_EQUAL(upw[2], expected_upw[2]); From 7e8ed40714851ca9343c7878df48a87f772d2151 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 1 Jun 2017 14:15:31 +0200 Subject: [PATCH 4/5] Use free function for multi-phase upwinding. --- opm/autodiff/BlackoilTransportModel.hpp | 55 +++++-------------------- 1 file changed, 10 insertions(+), 45 deletions(-) diff --git a/opm/autodiff/BlackoilTransportModel.hpp b/opm/autodiff/BlackoilTransportModel.hpp index d4d5bcc6d..41ddc56d6 100644 --- a/opm/autodiff/BlackoilTransportModel.hpp +++ b/opm/autodiff/BlackoilTransportModel.hpp @@ -26,6 +26,7 @@ #include #include #include +#include namespace Opm { @@ -401,58 +402,22 @@ namespace Opm { Eigen::Array multiPhaseUpwind(const std::vector& head_diff, const V& transmissibility) { - // Based on the paper "Upstream Differencing for Multiphase Flow in Reservoir Simulation", - // by Yann Brenier and Jérôme Jaffré, - // SIAM J. Numer. Anal., 28(3), 685–696. - // DOI:10.1137/0728036 - - // Using the data members: - // total_flux_ - // sd_.rq[].mob - - // Notation based on paper cited above. + assert(numPhases() == 3); const int num_connections = head_diff[0].size(); Eigen::Array upwind(num_connections, numPhases()); - using ValueAndIndex = std::pair; - const int num_phases = numPhases(); - std::vector g(num_phases); - std::vector theta(num_phases); for (int conn = 0; conn < num_connections; ++conn) { const double q = total_flux_[conn]; const double t = transmissibility[conn]; const int a = ops_.connection_cells(conn, 0); // first cell of connection const int b = ops_.connection_cells(conn, 1); // second cell of connection - - // Get and sort the g values (also called "weights" in the paper) for this connection. - for (int phase_idx = 0; phase_idx < num_phases; ++phase_idx) { - g[phase_idx] = ValueAndIndex(head_diff[phase_idx].value()[conn], phase_idx); - } - std::sort(g.begin(), g.end()); - - // Compute theta and r. - // Paper notation: subscript l -> ell (for read/searchability) - // Note that since we index phases from 0, r is one less than in the paper. - int r = -1; - for (int ell = 0; ell < num_phases; ++ell) { - theta[ell] = q; - for (int j = 0; j < num_phases; ++j) { - if (j < ell) { - theta[ell] += t * (g[ell].first - g[j].first) * sd_.rq[g[j].second].mob.value()[b]; - } - if (j > ell) { - theta[ell] += t * (g[ell].first - g[j].first) * sd_.rq[g[j].second].mob.value()[a]; - } - } - if (theta[ell] <= 0.0) { - r = ell; - } else { - break; // r is correct, no need to continue - } - } - - for (int ell = 0; ell < num_phases; ++ell) { - const int phase_idx = g[ell].second; - upwind(conn, phase_idx) = ell > r ? 1.0 : -1.0; + auto up = connectionMultiPhaseUpwind( + {{ head_diff[0].value()[conn], head_diff[1].value()[conn], head_diff[2].value()[conn] }}, + {{ sd_.rq[0].mob.value()[a], sd_.rq[1].mob.value()[a], sd_.rq[2].mob.value()[a]}}, + {{ sd_.rq[0].mob.value()[b], sd_.rq[1].mob.value()[b], sd_.rq[2].mob.value()[b]}}, + t, + q); + for (int ii = 0; ii < numPhases(); ++ii) { + upwind(conn, ii) = up[ii]; } } return upwind; From f5795322fa0cbbae0caa615c7f0bc8a3e114314d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 2 Jun 2017 08:52:16 +0200 Subject: [PATCH 5/5] Add missing include directive for std::sort(). --- opm/autodiff/multiPhaseUpwind.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/multiPhaseUpwind.cpp b/opm/autodiff/multiPhaseUpwind.cpp index 934208f50..043fe689e 100644 --- a/opm/autodiff/multiPhaseUpwind.cpp +++ b/opm/autodiff/multiPhaseUpwind.cpp @@ -20,6 +20,7 @@ #include +#include #include