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] 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