diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index a7b7a1b1a..cba690427 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -35,6 +35,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/BlackoilPropsAdFromDeck.cpp
+ opm/autodiff/WellDensitySegmented.cpp
)
# originally generated with the command:
@@ -45,6 +46,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_span.cpp
tests/test_syntax.cpp
tests/test_scalar_mult.cpp
+ tests/test_welldensitysegmented.cpp
)
list (APPEND TEST_DATA_FILES
@@ -86,4 +88,5 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorIncompTwophaseAd.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp
+ opm/autodiff/WellDensitySegmented.hpp
)
diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp
new file mode 100644
index 000000000..5ba38e87b
--- /dev/null
+++ b/opm/autodiff/WellDensitySegmented.cpp
@@ -0,0 +1,172 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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
+#include
+#include
+#include
+#include
+
+
+std::vector Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
+ const WellState& wstate,
+ const PhaseUsage& phase_usage,
+ const std::vector& b_perf,
+ const std::vector& rsmax_perf,
+ const std::vector& rvmax_perf,
+ const std::vector& z_perf,
+ const std::vector& surf_dens,
+ const double gravity)
+{
+ // Verify that we have consistent input.
+ const int np = wells.number_of_phases;
+ const int nw = wells.number_of_wells;
+ const int nperf = wells.well_connpos[nw];
+ if (wells.number_of_phases != phase_usage.num_phases) {
+ OPM_THROW(std::logic_error, "Inconsistent input: wells vs. phase_usage.");
+ }
+ if (surf_dens.size() != size_t(wells.number_of_phases)) {
+ OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage.");
+ }
+ if (nperf*np != int(wstate.perfRates().size())) {
+ OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
+ }
+ if (nperf*np != int(b_perf.size())) {
+ OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf.");
+ }
+ if (nperf != int(z_perf.size())) {
+ OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf.");
+ }
+ if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) {
+ // Need both oil and gas phases.
+ if (!phase_usage.phase_used[BlackoilPhases::Liquid]) {
+ OPM_THROW(std::logic_error, "Oil phase inactive, but non-empty rsmax_perf or rvmax_perf.");
+ }
+ if (!phase_usage.phase_used[BlackoilPhases::Vapour]) {
+ OPM_THROW(std::logic_error, "Gas phase inactive, but non-empty rsmax_perf or rvmax_perf.");
+ }
+ }
+
+ // Algorithm:
+
+ // We'll assume the perforations are given in order from top to
+ // bottom for each well. By top and bottom we do not necessarily
+ // mean in a geometric sense (depth), but in a topological sense:
+ // the 'top' perforation is nearest to the surface topologically.
+ // Our goal is to compute a pressure delta for each perforation.
+
+ // 1. Compute the flow (in surface volume units for each
+ // component) exiting up the wellbore from each perforation,
+ // taking into account flow from lower in the well, and
+ // in/out-flow at each perforation.
+ std::vector q_out_perf(nperf*np);
+ for (int w = 0; w < nw; ++w) {
+ // Iterate over well perforations from bottom to top.
+ for (int perf = wells.well_connpos[w+1] - 1; perf >= wells.well_connpos[w]; --perf) {
+ for (int phase = 0; phase < np; ++phase) {
+ if (perf == wells.well_connpos[w+1] - 1) {
+ // This is the bottom perforation. No flow from below.
+ q_out_perf[perf*np + phase] = 0.0;
+ } else {
+ // Set equal to flow from below.
+ q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
+ }
+ // Subtract outflow through perforation.
+ q_out_perf[perf*np + phase] -= wstate.perfRates()[perf*np + phase];
+ }
+ }
+ }
+
+ // 2. Compute the component mix at each perforation as the
+ // absolute values of the surface rates divided by their sum.
+ // Then compute volume ratios (formation factors) for each perforation.
+ // Finally compute densities for the segments associated with each perforation.
+ const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
+ const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
+ std::vector mix(np);
+ std::vector x(np);
+ std::vector dens(nperf);
+ for (int w = 0; w < nw; ++w) {
+ for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
+ // Find component mix.
+ const double tot_surf_rate = std::accumulate(q_out_perf.begin() + np*perf,
+ q_out_perf.begin() + np*(perf+1), 0.0);
+ if (tot_surf_rate != 0.0) {
+ for (int phase = 0; phase < np; ++phase) {
+ mix[phase] = std::fabs(q_out_perf[perf*np + phase]/tot_surf_rate);
+ }
+ } else {
+ // No flow => use well specified fractions for mix.
+ std::copy(wells.comp_frac + w*np, wells.comp_frac + (w+1)*np, mix.begin());
+ }
+ // Compute volume ratio.
+ x = mix;
+ double rs = 0.0;
+ double rv = 0.0;
+ if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
+ rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
+ }
+ if (!rvmax_perf.empty() && mix[gaspos] > 0.0) {
+ rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
+ }
+ if (rs != 0.0) {
+ // Subtract gas in oil from gas mixture
+ x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
+ }
+ if (rv != 0.0) {
+ // Subtract oil in gas from oil mixture
+ x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
+ }
+ double volrat = 0.0;
+ for (int phase = 0; phase < np; ++phase) {
+ volrat += x[phase] / b_perf[perf*np + phase];
+ }
+ // Compute segment density.
+ dens[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
+ }
+ }
+
+ // 3. Compute pressure differences between perforations.
+ // dp_perf will contain the pressure difference between a
+ // perforation and the one above it, except for the first
+ // perforation for each well, for which it will be the
+ // difference to the reference (bhp) depth.
+ std::vector dp_perf(nperf);
+ for (int w = 0; w < nw; ++w) {
+ for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) {
+ const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1];
+ const double dz = z_perf[perf] - z_above;
+ dp_perf[perf] = dz * dens[perf] * gravity;
+ }
+ }
+
+ // 4. Compute pressure differences to the reference point (bhp) by
+ // accumulating the already computed adjacent pressure
+ // differences, storing the result in dp_perf.
+ // This accumulation must be done per well.
+ for (int w = 0; w < nw; ++w) {
+ const auto beg = dp_perf.begin() + wells.well_connpos[w];
+ const auto end = dp_perf.begin() + wells.well_connpos[w + 1];
+ std::partial_sum(beg, end, beg);
+ }
+
+ return dp_perf;
+}
diff --git a/opm/autodiff/WellDensitySegmented.hpp b/opm/autodiff/WellDensitySegmented.hpp
new file mode 100644
index 000000000..228262311
--- /dev/null
+++ b/opm/autodiff/WellDensitySegmented.hpp
@@ -0,0 +1,66 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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_WELLDENSITYSEGMENTED_HEADER_INCLUDED
+#define OPM_WELLDENSITYSEGMENTED_HEADER_INCLUDED
+
+#include
+
+struct Wells;
+
+namespace Opm
+{
+
+ class WellState;
+ class PhaseUsage;
+
+
+ /// A class giving a well model, by which we mean a way to compute
+ /// the pressure deltas of each perforation and the bottom-hole
+ /// pressure. This class contains an explicit model, that uses a
+ /// different density for each well segment, that is between each
+ /// pair of perforations.
+ class WellDensitySegmented
+ {
+ public:
+ /// Compute pressure deltas.
+ /// Notation: N = number of perforations, P = number of phases.
+ /// \param[in] wells struct with static well info
+ /// \param[in] wstate dynamic well solution information, only perfRates() is used
+ /// \param[in] phase_usage specifies which phases are active and not
+ /// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation
+ /// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N
+ /// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N
+ /// \param[in] z_perf depth values for each perforation, size N
+ /// \param[in] surf_dens surface densities for active components, size P
+ /// \param[in] gravity gravity acceleration constant
+ static std::vector computeConnectionPressureDelta(const Wells& wells,
+ const WellState& wstate,
+ const PhaseUsage& phase_usage,
+ const std::vector& b_perf,
+ const std::vector& rsmax_perf,
+ const std::vector& rvmax_perf,
+ const std::vector& z_perf,
+ const std::vector& surf_dens,
+ const double gravity);
+ };
+
+} // namespace Opm
+
+#endif // OPM_WELLDENSITYSEGMENTED_HEADER_INCLUDED
diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp
new file mode 100644
index 000000000..a845c0fde
--- /dev/null
+++ b/tests/test_welldensitysegmented.cpp
@@ -0,0 +1,101 @@
+/*
+ Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+ 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 BOOST_TEST_MODULE WellDensitySegmentedTest
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#include
+#include
+
+using namespace Opm;
+
+BOOST_AUTO_TEST_CASE(TestPressureDeltas)
+{
+ // Simple water injector.
+ const int np = 3;
+ const int nperf = 10;
+ const double ref_depth = 0.0;
+ const double comp_frac_w[np] = { 1.0, 0.0, 0.0 };
+ const double comp_frac_o[np] = { 0.0, 1.0, 0.0 };
+ const int cells[nperf/2] = { 0, 1, 2, 3, 4 };
+ const double WI[nperf/2] = { 1.0, 1.0, 1.0, 1.0, 1.0 };
+ std::shared_ptr wells(create_wells(np, 2, nperf), destroy_wells);
+ BOOST_REQUIRE(wells);
+ int ok = add_well(INJECTOR, ref_depth, nperf/2, comp_frac_w, cells, WI, "INJ", wells.get());
+ BOOST_REQUIRE(ok);
+ ok = add_well(PRODUCER, ref_depth, nperf/2, comp_frac_o, cells, WI, "PROD", wells.get());
+ BOOST_REQUIRE(ok);
+ std::vector rates = { 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0 };
+ WellState wellstate;
+ wellstate.perfRates() = rates;
+ PhaseUsage pu;
+ pu.num_phases = 3;
+ pu.phase_used[0] = true;
+ pu.phase_used[1] = true;
+ pu.phase_used[2] = true;
+ pu.phase_pos[0] = 0;
+ pu.phase_pos[1] = 1;
+ pu.phase_pos[2] = 2;
+ const std::vector b_perf = { 2.0, 3.0, 100,
+ 2.1, 3.3, 110,
+ 2.2, 3.6, 120,
+ 2.3, 4.0, 130,
+ 2.4, 4.5, 140,
+ 2.0, 3.0, 100,
+ 2.1, 3.3, 110,
+ 2.2, 3.6, 120,
+ 2.3, 4.0, 130,
+ 2.4, 4.5, 140 };
+ const std::vector rsmax_perf = { 50, 50, 50, 50, 50, 50, 50, 50, 50, 50 };
+ const std::vector rvmax_perf = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ const std::vector z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 };
+ const std::vector surf_dens = { 1000.0, 800.0, 10.0 };
+ const double gravity = Opm::unit::gravity;
+ const std::vector dp = WellDensitySegmented::computeConnectionPressureDelta(*wells, wellstate, pu, b_perf, rsmax_perf, rvmax_perf, z_perf, surf_dens, gravity);
+ const std::vector answer = { 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity,
+ 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity };
+ BOOST_REQUIRE_EQUAL(dp.size(), answer.size());
+ // for (auto p : dp) { std::cout << p << std::endl; }
+ for (size_t i = 0; i < dp.size(); ++i) {
+ BOOST_CHECK_CLOSE(dp[i], answer[i], 1e-8);
+ }
+}