From 491b532a6a14f939557a7dc8ca1557fd8f368830 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Tue, 18 Aug 2020 09:13:54 +0200
Subject: [PATCH] Add trivial saltvd table for cases without saltvd given in
 the deck

---
 ebos/equil/initstateequil.hh | 30 +++++++++++++++++++++---------
 tests/test_equil.cc          | 24 ++++++++++++++++++++++++
 2 files changed, 45 insertions(+), 9 deletions(-)

diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh
index a4a731f41..8cabb11db 100644
--- a/ebos/equil/initstateequil.hh
+++ b/ebos/equil/initstateequil.hh
@@ -178,7 +178,8 @@ private:
     density(const double depth,
             const double press) const
     {
-        double saltConcentration = saltVdTable_.eval(depth);
+        // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate
+        double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true);
         double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration);
         rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
         return rho;
@@ -1736,17 +1737,28 @@ private:
     template <class RMap>
     void updateInitialSaltConcentration_(const Opm::EclipseState& eclState, const RMap& reg, const Grid& grid)
     {
+        const int numEquilReg = rsFunc_.size();
+        saltVdTable_.resize(numEquilReg);
         const auto& tables = eclState.getTableManager();
         const Opm::TableContainer& saltvdTables = tables.getSaltvdTables();
-        saltVdTable_.resize(saltvdTables.size());
-        for (size_t i = 0; i < saltvdTables.size(); ++i) {
-            const Opm::SaltvdTable& saltvdTable = saltvdTables.getTable<Opm::SaltvdTable>(i);
-            saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
 
-            const auto& cells = reg.cells(i);
-            for (const auto& cell : cells) {
-                const double depth = UgGridHelpers::cellCenterDepth(grid, cell);
-                this->saltConcentration_[cell] = saltVdTable_[i].eval(depth);
+        // If no saltvd table is given, we create a trivial table for the density calculations
+        if (saltvdTables.empty()) {
+            std::vector<double> x = {0.0,1.0};
+            std::vector<double> y = {0.0,0.0};
+            for (size_t i = 0; i < numEquilReg; ++i) {
+                saltVdTable_[i].setXYContainers(x,y);
+            }
+        } else {
+            for (size_t i = 0; i < saltvdTables.size(); ++i) {
+                const Opm::SaltvdTable& saltvdTable = saltvdTables.getTable<Opm::SaltvdTable>(i);
+                saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
+
+                const auto& cells = reg.cells(i);
+                for (const auto& cell : cells) {
+                    const double depth = UgGridHelpers::cellCenterDepth(grid, cell);
+                    this->saltConcentration_[cell] = saltVdTable_[i].eval(depth);
+                }
             }
         }
     }
diff --git a/tests/test_equil.cc b/tests/test_equil.cc
index a278ff64c..071ba4646 100644
--- a/tests/test_equil.cc
+++ b/tests/test_equil.cc
@@ -182,6 +182,11 @@ void test_PhasePressure()
     using TypeTag     = TTAG(TestEquilTypeTag);
     using FluidSystem = GET_PROP_TYPE(TypeTag, FluidSystem);
 
+    using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
+    std::vector<double> x = {0.0,100.0};
+    std::vector<double> y = {0.0,0.0};
+    TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
+
     auto simulator = initSimulator<TypeTag>("equil_base.DATA");
     initDefaultFluidSystem<TypeTag>();
 
@@ -189,6 +194,7 @@ void test_PhasePressure()
         record,
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0
     };
 
@@ -235,26 +241,35 @@ void test_CellSubset()
     const Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
                                   mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
 
+    using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
+    std::vector<double> x = {0.0,100.0};
+    std::vector<double> y = {0.0,0.0};
+    TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
+
     const Opm::EQUIL::EquilReg region[] =
     {
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
     };
 
@@ -330,26 +345,35 @@ void test_RegMapping()
     auto simulator = initSimulator<TypeTag>("equil_base.DATA");
     initDefaultFluidSystem<TypeTag>();
 
+    using TabulatedFunction = typename Opm::Tabulated1DFunction<double>;
+    std::vector<double> x = {0.0,100.0};
+    std::vector<double> y = {0.0,0.0};
+    TabulatedFunction trivialSaltVdTable = TabulatedFunction(2, x, y);
+
     const Opm::EQUIL::EquilReg region[] =
     {
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
+        trivialSaltVdTable,
         0)
         };