Merge pull request #611 from totto82/TLpressure

Implement pressure effects in the Todd-Longstaff mixing parameter

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]  po              Array of n oil pressure values
+    /// \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/ParseContext.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::ParseContext parseContext;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseContext);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
+        const int global_ind = 0;
+        Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 1, &global_ind);
+}
+
+BOOST_AUTO_TEST_CASE(SolventData)
+{
+        Opm::ParseContext parseContext;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseContext);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
+        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::ParseContext parseContext;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData + pmiscData, parseContext);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
+        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::ParseContext parseContext;
+        Opm::ParserPtr parser(new Opm::Parser());
+        Opm::DeckPtr deck =  parser->parseString(deckData + solventData + tlpmixpaData, parseContext);
+        Opm::EclipseStateConstPtr eclState;
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
+        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::ParseContext parseContext;
+    Opm::ParserPtr parser(new Opm::Parser());
+    // no pmisc data and default tlpmixdata i.e it should throw
+    Opm::DeckPtr deck =  parser->parseString(deckData + solventData, parseContext);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseContext));
+    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::ParseContext parseContext;
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::DeckPtr deck =  parser->parseString(deckData + solventData + pmiscData + tlpmixpaDataDefault, parseContext);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseContext));
+    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::ParseContext parseContext;
+    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, parseContext);
+    Opm::EclipseStateConstPtr eclState;
+    eclState.reset(new Opm::EclipseState(deck , parseContext));
+    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);
+}
+