Implement pressure effects in the Todd-Longstaff mixing parameter

The Todd-Longstaff model is extended to incorporate pressure effects
The solvent viscosity is then caculated as

mu_eff = mu_s^(1-\alpha * \omega) * mu_mix^(\alpha * \omega)

where \omega accounts for the porous media effects and \alpha =
\alpha(pressure) accounts for the miscibility of the solvent and oil
when contacted.
The \alpha values can be given using the TLPMIXPA keyword

If no entries are given to TLPMIXPA the table specified using PMISC will
be used as default.
IF TLPMIXPA does not appear in the grid \alpha = 1 and the pressure
effect is neglected.
This is tested in test_solventprops_ad.cpp

CMakeLists_files.cmake

@@ -79,6 +79,7 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_welldensitysegmented.cpp
 	tests/test_vfpproperties.cpp
 	tests/test_singlecellsolves.cpp
+	tests/test_solventprops_ad.cpp
 	)
 
 list (APPEND TEST_DATA_FILES

opm/autodiff/BlackoilSolventModel.hpp

@@ -233,7 +233,7 @@ namespace Opm {
         void computeEffectiveProperties(const SolutionState&  state);
 
         // compute density and viscosity using the ToddLongstaff mixing model
-        void computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu);
+        void computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const ADB po, const Opm::PhaseUsage pu);
 
         // compute phase pressures.
         std::vector<ADB>

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -828,7 +828,7 @@ namespace Opm {
         effective_saturations[solvent_pos_] = ss - sgcwmis;
 
         // Compute effective viscosities and densities
-        computeToddLongstaffMixing(viscosity, density, effective_saturations, pu);
+        computeToddLongstaffMixing(viscosity, density, effective_saturations, po, pu);
 
         // compute the volume factors from the densities
         const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ],  cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs);
@@ -855,7 +855,7 @@ namespace Opm {
 
     template <class Grid>
     void
-    BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu)
+    BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const ADB po, const Opm::PhaseUsage pu)
     {
         const int  nc = cells_.size();
         const V ones = V::Constant(nc, 1.0);
@@ -897,7 +897,9 @@ namespace Opm {
         const ADB mu_m = zero_selectorSn.select(mu_s + mu_o + mu_g, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow *  mu_g_pow)
                                                        + ( (ss_eff / sn_eff) * mu_o_pow *  mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
         // Mixing parameter for viscosity
-        const V mix_param_mu = solvent_props_.mixingParameterViscosity(cells_);
+        // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterViscosity the effect of the porous media.
+        // The pressureMixingParameter is not implemented in ecl100.
+        const ADB mix_param_mu = solvent_props_.mixingParameterViscosity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
 
         // Update viscosities
         viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,ones - mix_param_mu) * pow(mu_mos, mix_param_mu);
@@ -910,7 +912,7 @@ namespace Opm {
         ADB& rho_s = density[solvent_pos_];
 
         // mixing parameter for density
-        const V mix_param_rho = solvent_props_.mixingParameterDensity(cells_);
+        const ADB mix_param_rho = solvent_props_.mixingParameterDensity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
 
         // compute effective viscosities for density calculations. These have to
         // be recomputed as a different mixing parameter may be used.

opm/autodiff/SolventPropsAdFromDeck.cpp

@@ -27,6 +27,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/PvdsTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
+
 
 namespace Opm
 {
@@ -274,6 +276,35 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
                 }
             }
 
+            if (deck->hasKeyword("TLPMIXPA")) {
+                const TableContainer& tlpmixparTables = tables->getTlpmixpaTables();
+                if (!tlpmixparTables.empty()) {
+
+                    int numRegions = tlpmixparTables.size();
+                    // resize the attributes of the object
+                    tlpmix_param_.resize(numRegions);
+                    for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+                        const Opm::TlpmixpaTable& tlpmixparTable = tlpmixparTables.getTable<TlpmixpaTable>(regionIdx);
+
+                        // Copy data
+                        const auto& po = tlpmixparTable.getOilPhasePressureColumn();
+                        const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn();
+
+                        tlpmix_param_[regionIdx] = NonuniformTableLinear<double>(po, tlpmixpa);
+
+                    }
+                } else {
+                    // if empty keyword. Try to use the pmisc table as default.
+                    if (pmisc_.size() > 0) {
+                        tlpmix_param_ = pmisc_;
+                    } else {
+                        OPM_THROW(std::invalid_argument, "If the pressure dependent TL values in TLPMIXPA is defaulted (no entries), then the PMISC tables must be specified.");
+                    }
+                }
+            }
+
+
+
         }
     }
 
@@ -451,6 +482,15 @@ V SolventPropsAdFromDeck::mixingParameterDensity(const Cells& cells) const {
     return V::Zero(n);
 }
 
+ADB SolventPropsAdFromDeck::pressureMixingParameter(const ADB& po,
+                                                    const Cells& cells) const {
+    if (tlpmix_param_.size() > 0) {
+        return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, tlpmix_param_);
+    }
+    // return ones if not specified i.e. no pressure effects.
+    return ADB::constant(V::Constant(po.size(), 1.0));
+}
+
 void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword,
                                                Opm::EclipseStateConstPtr eclState,
                                                size_t numCompressed,

opm/autodiff/SolventPropsAdFromDeck.hpp

@@ -143,6 +143,13 @@ public:
     /// return                      Array of n mixing paramters for density calculation
     V mixingParameterDensity(const Cells& cells) const;
 
+    /// Todd-Longstaff pressure dependent mixing parameter
+    /// \param[in]  So              Array of n oil fraction values. Soil / Sn values, where Sn = Sgas + Ssolvent + Soil.
+    /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
+    /// return                      Array of n pressure dependent mixing paramters
+    ADB pressureMixingParameter(const ADB& po,
+                              const Cells& cells) const;
+
 
 private:
     /// Makes ADB from table values
@@ -185,6 +192,7 @@ private:
     std::vector<NonuniformTableLinear<double> > pmisc_;
     std::vector<NonuniformTableLinear<double> > sorwmis_;
     std::vector<NonuniformTableLinear<double> > sgcwmis_;
+    std::vector<NonuniformTableLinear<double> > tlpmix_param_;
     std::vector<double> mix_param_viscosity_;
     std::vector<double> mix_param_density_;
 };

tests/test_solventprops_ad.cpp

@@ -0,0 +1,247 @@
+/*
+  Copyright 2016 IRIS 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 <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define BOOST_TEST_MODULE SolventPropertiesTest
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
+
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <fstream>
+#include <iostream>
+
+
+const std::string deckData = "\n\
+RUNSPEC \n\
+\n\
+SOLVENT \n\
+\n\
+MISCIBLE\n\
+1  3 /\n\
+\n\
+DIMENS \n\
+1 1 1 \n\
+/\n\
+TABDIMS\n\
+/\n\
+GRID \n\
+\n\
+DXV \n\
+1 \n\
+/\n\
+DYV \n\
+1 \n\
+/\n\
+DZV \n\
+1 \n\
+/\n";
+
+const std::string solventData = "\n\
+SDENSITY \n\
+0.1 / \n\
+PVDS \n\
+1 1 0.1 / \n\
+SSFN \n\
+0.0 0.0 0.0 \n\
+1.0 1.0 1.0 \n\
+/ \n\
+MISC \n\
+0.0 0.0 \n\
+1.0 1.0 \n\
+/ \n\
+SOF2 \n\
+0 0  \n\
+0.88 1 / \n";
+
+BOOST_AUTO_TEST_CASE(Construction)
+{
+        Opm::ParseMode parseMode;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseMode);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        const int global_ind = 0;
+        Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 1, &global_ind);
+}
+
+BOOST_AUTO_TEST_CASE(SolventData)
+{
+        Opm::ParseMode parseMode;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseMode);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        const int global_ind = 0;
+        Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 1, &global_ind);
+
+        const Opm::SolventPropsAdFromDeck::Cells cells(1, 0);
+        typedef Opm::SolventPropsAdFromDeck::V V;
+        V rho = solventprops.solventSurfaceDensity(cells);
+        BOOST_REQUIRE_EQUAL(rho.size(), cells.size());
+        BOOST_CHECK_EQUAL(rho[0], 0.1);
+}
+
+const std::string pmiscData = "\n\
+PMISC\n\
+100 0.0 \n\
+200 0.0 \n\
+500 1.0 \n\
+1000 1.0 /\n\
+\n";
+
+BOOST_AUTO_TEST_CASE(PMISC)
+{
+        Opm::ParseMode parseMode;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData + pmiscData, parseMode);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
+        typedef Opm::SolventPropsAdFromDeck::V V;
+        const int* global_ind = new int[3] {0 , 1 , 2};
+        Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
+        V po(3);
+        po << 150,250,550;
+        po = po * Opm::unit::barsa;
+        BOOST_REQUIRE_EQUAL(po.size(), cells.size());
+        V pmisc = solventprops.pressureMiscibilityFunction(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
+        BOOST_REQUIRE_EQUAL(pmisc.size(), cells.size());
+        BOOST_CHECK_EQUAL(pmisc[0], 0.0);
+        const double tol = 1e-12;
+        const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
+        BOOST_CHECK_CLOSE(pmisc[1], value, tol);
+        BOOST_CHECK_EQUAL(pmisc[2], 1.0);
+}
+
+const std::string tlpmixpaData = "\n\
+TLPMIXPA\n\
+100 0.0 \n\
+200 0.0 \n\
+500 1.0 \n\
+1000 1.0 /\n\
+\n";
+
+
+BOOST_AUTO_TEST_CASE(TLPMIXPA)
+{
+        Opm::ParseMode parseMode;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData + tlpmixpaData, parseMode);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
+        typedef Opm::SolventPropsAdFromDeck::V V;
+        const int* global_ind = new int[3] {0 , 1 , 2};
+        Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
+        V po(3);
+        po << 150,250,550;
+        po = po * Opm::unit::barsa;
+        BOOST_REQUIRE_EQUAL(po.size(), cells.size());
+        V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
+        BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
+        BOOST_CHECK_EQUAL(tlpmixpa[0], 0.0);
+        const double tol = 1e-12;
+        const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
+        BOOST_CHECK_CLOSE(tlpmixpa[1], value, tol);
+        BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
+}
+
+BOOST_AUTO_TEST_CASE(TLPMIXPA_NOT_SPECIFIED)
+{
+    Opm::ParseMode parseMode;
+    Opm::ParserPtr parser(new Opm::Parser());
+    // no pmisc data and default tlpmixdata i.e it should throw
+    Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseMode);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseMode));
+    const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
+    const int* global_ind = new int[3] {0 , 1 , 2};
+    Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
+    typedef Opm::SolventPropsAdFromDeck::V V;
+    V po(3);
+    po << 150,250,550;
+    po = po * Opm::unit::barsa;
+    BOOST_REQUIRE_EQUAL(po.size(), cells.size());
+    V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
+    BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
+    // if not specified tlpmixpa is 1.0 for all cells.
+    BOOST_CHECK_EQUAL(tlpmixpa[0], 1.0);
+    BOOST_CHECK_EQUAL(tlpmixpa[1], 1.0);
+    BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
+
+}
+
+const std::string tlpmixpaDataDefault = "\n\
+TLPMIXPA\n\
+/\n\
+\n";
+
+BOOST_AUTO_TEST_CASE(TLPMIXPA_DEFAULT)
+{
+    Opm::ParseMode parseMode;
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::DeckPtr deck =  parser->parseString(deckData + solventData + pmiscData + tlpmixpaDataDefault, parseMode);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseMode));
+    const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
+    typedef Opm::SolventPropsAdFromDeck::V V;
+    const int* global_ind = new int[3] {0 , 1 , 2};
+    Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
+    V po(3);
+    po << 150,250,550;
+    po = po * Opm::unit::barsa;
+    BOOST_REQUIRE_EQUAL(po.size(), cells.size());
+    V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
+    BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
+    BOOST_CHECK_EQUAL(tlpmixpa[0], 0.0);
+    const double tol = 1e-12;
+    const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
+    BOOST_CHECK_CLOSE(tlpmixpa[1], value, tol);
+    BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
+}
+
+BOOST_AUTO_TEST_CASE(TLPMIXPA_DEFAULT_NOPMISC)
+{
+    Opm::ParseMode parseMode;
+    Opm::ParserPtr parser(new Opm::Parser());
+    // no pmisc data and default tlpmixdata i.e it should throw
+    Opm::DeckPtr deck =  parser->parseString(deckData + solventData + tlpmixpaDataDefault, parseMode);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseMode));
+    const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
+    const int* global_ind = new int[3] {0 , 1 , 2};
+    BOOST_CHECK_THROW(Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind), std::invalid_argument);
+}
+