diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index fcabf8bb6..d42d26037 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 @@ -90,6 +91,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 @@ -186,6 +188,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/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; diff --git a/opm/autodiff/multiPhaseUpwind.cpp b/opm/autodiff/multiPhaseUpwind.cpp new file mode 100644 index 000000000..043fe689e --- /dev/null +++ b/opm/autodiff/multiPhaseUpwind.cpp @@ -0,0 +1,88 @@ +/* + 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 +#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 diff --git a/tests/test_multiphaseupwind.cpp b/tests/test_multiphaseupwind.cpp new file mode 100644 index 000000000..cac7fc806 --- /dev/null +++ b/tests/test_multiphaseupwind.cpp @@ -0,0 +1,108 @@ +/* + 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 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. + +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 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(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]); +} + + +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 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(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]); +} + + +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 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(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]); +}