Add support for ssfn

The SSFN keyword is used in the solvent model to modify the gas/solvent
rel.perm values.

opm/autodiff/ExtendedBlackoilModel_impl.hpp

@@ -379,11 +379,11 @@ namespace Opm {
                 ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
                 V ones = V::Constant(nc, 1.0);
 
-                // TODO: Add support for gas/solvent function (SSFN)
                 const ADB tr_mult = transMult(state.pressure);
                 const ADB mu = solvent_props_.muSolvent(phasePressure,cells_);               
-                rq_[solvent_pos_].mob = F_solvent * tr_mult * kr / mu;
+
-                rq_[actph].mob = (ones - F_solvent) * rq_[actph].mob;
+                rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu;
+                rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob;
 
                 const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b;
                 const ADB rhoavg_solvent = ops_.caver * rho_solvent;

opm/autodiff/SolventPropsAdFromDeck.cpp

@@ -87,6 +87,34 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
     } else {
         OPM_THROW(std::runtime_error, "PVDS must be specified in SOLVENT runs\n");
     }
+
+    const auto& ssfnTables = eclState->getSsfnTables();
+    // relative permeabilty multiplier
+    if (!ssfnTables.empty()) {
+
+        int numRegions = pvdsTables.size();
+
+        if(numRegions > 1) {
+            OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN");
+        }
+        // resize the attributes of the object
+        krg_.resize(numRegions);
+        krs_.resize(numRegions);
+        for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+            const Opm::SsfnTable& ssfnTable = ssfnTables[regionIdx];
+
+            // Copy data
+            const std::vector<double>& solventFraction = ssfnTable.getSolventFractionColumn();
+            const std::vector<double>& krg = ssfnTable.getGasRelPermMultiplierColumn();
+            const std::vector<double>& krs = ssfnTable.getSolventRelPermMultiplierColumn();
+
+            krg_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krg);
+            krs_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krs);
+        }
+
+    } else {
+        OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n");
+    }
 }
 
 ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
@@ -139,6 +167,52 @@ ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg,
     return ADB::function(std::move(b), std::move(jacs));
 }
 
+ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction,
+                                 const Cells& cells) const
+{
+    const int n = cells.size();
+    assert(solventFraction.value().size() == n);
+    V krg(n);
+    V dkrgdsf(n);
+    for (int i = 0; i < n; ++i) {
+        const double& solventFraction_i = solventFraction.value()[i];
+        int regionIdx = 0; // TODO add mapping from cells to sat function table
+        krg[i] = krg_[regionIdx](solventFraction_i);
+        dkrgdsf[i] = krg_[regionIdx].derivative(solventFraction_i);
+    }
+
+    ADB::M dkrgdsf_diag = spdiag(dkrgdsf);
+    const int num_blocks = solventFraction.numBlocks();
+    std::vector<ADB::M> jacs(num_blocks);
+    for (int block = 0; block < num_blocks; ++block) {
+        fastSparseProduct(dkrgdsf_diag, solventFraction.derivative()[block], jacs[block]);
+    }
+    return ADB::function(std::move(krg), std::move(jacs));
+}
+
+ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction,
+                                 const Cells& cells) const
+{
+    const int n = cells.size();
+    assert(solventFraction.value().size() == n);
+    V krs(n);
+    V dkrsdsf(n);
+    for (int i = 0; i < n; ++i) {
+        const double& solventFraction_i = solventFraction.value()[i];
+        int regionIdx = 0; // TODO add mapping from cells to sat function table
+        krs[i] = krs_[regionIdx](solventFraction_i);
+        dkrsdsf[i] = krs_[regionIdx].derivative(solventFraction_i);
+    }
+
+    ADB::M dkrsdsf_diag = spdiag(dkrsdsf);
+    const int num_blocks = solventFraction.numBlocks();
+    std::vector<ADB::M> jacs(num_blocks);
+    for (int block = 0; block < num_blocks; ++block) {
+        fastSparseProduct(dkrsdsf_diag, solventFraction.derivative()[block], jacs[block]);
+    }
+    return ADB::function(std::move(krs), std::move(jacs));
+}
+
 V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
     const int n = cells.size();
     V density(n);

opm/autodiff/SolventPropsAdFromDeck.hpp

@@ -65,6 +65,20 @@ public:
     ADB muSolvent(const ADB& pg,
               const Cells& cells) const;
 
+    /// Gas relPerm multipliers
+    /// \param[in]  solventFraction Array of n gas pressure values.
+    /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
+    /// \return                     Array of n gas relPerm multiplier values.
+    ADB gasRelPermMultiplier(const ADB& solventFraction,
+              const Cells& cells) const;
+
+    /// Solvent relPerm multipliers
+    /// \param[in]  solventFraction Array of n gas pressure values.
+    /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
+    /// \return                     Array of n solvent relPerm multiplier values.
+    ADB solventRelPermMultiplier(const ADB& solventFraction,
+              const Cells& cells) const;
+
     V solventSurfaceDensity(const Cells& cells) const;
 
 private:
@@ -74,6 +88,8 @@ private:
     std::vector<NonuniformTableLinear<double> > viscosity_;
     std::vector<NonuniformTableLinear<double> > inverseBmu_;
     std::vector<double> solvent_surface_densities_;
+    std::vector<NonuniformTableLinear<double> > krg_;
+    std::vector<NonuniformTableLinear<double> > krs_;
 };
 
 } // namespace OPM