OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1199 from atgeirr/multiphase-upwind-refactor

Browse Source 
Multiphase upwind refactor
This commit is contained in:
Atgeirr Flø Rasmussen
2017-06-02 11:13:14 +02:00
committed by GitHub
parent 4c70c831d4 f5795322fa
commit 37cbdbced2
5 changed files with 254 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

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

55
opm/autodiff/BlackoilTransportModel.hpp
View File

@@ -26,6 +26,7 @@
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
#include <opm/autodiff/multiPhaseUpwind.hpp>
namespace Opm {
@@ -401,58 +402,22 @@ namespace Opm {
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> multiPhaseUpwind(const std::vector<ADB>& 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<double, Eigen::Dynamic, Eigen::Dynamic> upwind(num_connections, numPhases());
using ValueAndIndex = std::pair<double, int>;
const int num_phases = numPhases();
std::vector<ValueAndIndex> g(num_phases);
std::vector<double> 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;

88
opm/autodiff/multiPhaseUpwind.cpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/multiPhaseUpwind.hpp>
#include <algorithm>
#include <utility>
namespace Opm
{
std::array<double, 3> connectionMultiPhaseUpwind(const std::array<double, 3>& head_diff,
const std::array<double, 3>& mob1,
const std::array<double, 3>& 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<double, int>;
std::array<ValueAndIndex, NumPhases> 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<double, NumPhases> 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<double, NumPhases> 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

45
opm/autodiff/multiPhaseUpwind.hpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTIPHASEUPWIND_HEADER_INCLUDED
#define OPM_MULTIPHASEUPWIND_HEADER_INCLUDED
#include <array>
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<double, 3> connectionMultiPhaseUpwind(const std::array<double, 3>& head_diff,
const std::array<double, 3>& mob1,
const std::array<double, 3>& mob2,
const double transmissibility,
const double flux);
} // namespace Opm
#endif // OPM_MULTIPHASEUPWIND_HEADER_INCLUDED

108
tests/test_multiphaseupwind.cpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <opm/autodiff/multiPhaseUpwind.hpp>
// 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<double, 3> gd = {{ 4.0, -1.0, -2.0 }};
const std::array<double, 3> mob1 = {{ 1.0, 1.0, 1.0 }};
const std::array<double, 3> mob2 = {{ 1.0, 1.0, 1.0 }};
const double transmissibility = 1.0;
const double flux = 1.0;
const std::array<double, 3> expected_upw = {{ 1.0, -1.0, -1.0 }};
const std::array<double, 3> 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<double, 3> gd = {{ 4.0, -1.0, -2.0 }};
const std::array<double, 3> mob1 = {{ 1.0, 1.0, 1.0 }};
const std::array<double, 3> mob2 = {{ 1.0, 1.0, 1.0 }};
const double transmissibility = 1.0;
const double flux = 5.0;
const std::array<double, 3> expected_upw = {{ 1.0, 1.0, -1.0 }};
const std::array<double, 3> 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<double, 3> gd = {{ 4.0, -1.0, -2.0 }};
const std::array<double, 3> mob1 = {{ 1.0, 1.0, 1.0 }};
const std::array<double, 3> mob2 = {{ 1.0, 1.0, 1.0 }};
const double transmissibility = 1.0;
const double flux = 10.0;
const std::array<double, 3> expected_upw = {{ 1.0, 1.0, 1.0 }};
const std::array<double, 3> 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]);
}