Merge pull request #3759 from plgbrts/saltp

Enable salt precipitation/dissolution

CMakeLists.txt

@@ -388,7 +388,7 @@ add_dependencies(moduleVersion opmsimulators)
 set(FLOW_MODELS blackoil brine energy extbo foam gasoil gaswater
                 oilwater oilwater_brine gaswater_brine oilwater_polymer
                 oilwater_polymer_injectivity micp polymer solvent
-                gasoil_energy)
+                gasoil_energy brine_saltprecipitation)
 set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy)
 
 set(FLOW_TGTS)

ebos/eclequilinitializer.hh

@@ -153,9 +153,13 @@ public:
 
             }
 
-            // set salt concentration
+            // set the salt concentration
             if (enableBrine)
                 fluidState.setSaltConcentration(initialState.saltConcentration()[elemIdx]);
+                
+            //set the (solid) salt saturation
+            if (enableSaltPrecipitation)
+                fluidState.setSaltSaturation(initialState.saltSaturation()[elemIdx]);
         }
     }
 

ebos/eclgenericoutputblackoilmodule.cc

@@ -75,6 +75,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState,
                            bool enablePolymer,
                            bool enableFoam,
                            bool enableBrine,
+                           bool enableSaltPrecipitation,
                            bool enableExtbo,
                            bool enableMICP)
     : eclState_(eclState)
@@ -87,6 +88,7 @@ EclGenericOutputBlackoilModule(const EclipseState& eclState,
     , enablePolymer_(enablePolymer)
     , enableFoam_(enableFoam)
     , enableBrine_(enableBrine)
+    , enableSaltPrecipitation_(enableSaltPrecipitation)
     , enableExtbo_(enableExtbo)
     , enableMICP_(enableMICP)
 {
@@ -655,6 +657,8 @@ assignToSolution(data::Solution& sol)
         {"RV",       UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_SOLUTION,  rv_},
         {"RVSAT",    UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_AUXILIARY, oilVaporizationFactor_},
         {"SALT",     UnitSystem::measure::salinity,  data::TargetType::RESTART_SOLUTION,      cSalt_},
+        {"SALTP",    UnitSystem::measure::identity,  data::TargetType::RESTART_AUXILIARY,     pSalt_},
+        {"PERMFACT", UnitSystem::measure::identity,  data::TargetType::RESTART_AUXILIARY,     permFact_},
         {"SOMAX",    UnitSystem::measure::identity,  data::TargetType::RESTART_SOLUTION,      soMax_},
         {"SSOLVENT", UnitSystem::measure::identity,  data::TargetType::RESTART_SOLUTION,      sSol_},
         {"SS_X",     UnitSystem::measure::identity,  data::TargetType::RESTART_SOLUTION,      extboX_},
@@ -761,6 +765,10 @@ setRestart(const data::Solution& sol,
         cFoam_[elemIdx] = sol.data("FOAM")[globalDofIndex];
     if (!cSalt_.empty() && sol.has("SALT"))
         cSalt_[elemIdx] = sol.data("SALT")[globalDofIndex];
+    if (!pSalt_.empty() && sol.has("SALTP"))
+        pSalt_[elemIdx] = sol.data("SALTP")[globalDofIndex];
+    if (!permFact_.empty() && sol.has("PERMFACT"))
+        permFact_[elemIdx] = sol.data("PERMFACT")[globalDofIndex];
     if (!soMax_.empty() && sol.has("SOMAX"))
         soMax_[elemIdx] = sol.data("SOMAX")[globalDofIndex];
     if (!pcSwMdcOw_.empty() &&sol.has("PCSWM_OW"))
@@ -954,6 +962,10 @@ doAllocBuffers(unsigned bufferSize,
         cFoam_.resize(bufferSize, 0.0);
     if (enableBrine_)
         cSalt_.resize(bufferSize, 0.0);
+    if (enableSaltPrecipitation_) {
+        pSalt_.resize(bufferSize, 0.0);
+        permFact_.resize(bufferSize, 0.0);
+    }
     if (enableExtbo_) {
         extboX_.resize(bufferSize, 0.0);
         extboY_.resize(bufferSize, 0.0);

ebos/eclgenericoutputblackoilmodule.hh

@@ -128,6 +128,22 @@ public:
         return 0;
     }
 
+    Scalar getSaltSaturation(unsigned elemIdx) const
+    {
+        if (pSalt_.size() > elemIdx)
+            return pSalt_[elemIdx];
+
+        return 0;
+    }
+
+    Scalar getPermFactor(unsigned elemIdx) const
+    {
+        if (permFact_.size() > elemIdx)
+            return permFact_[elemIdx];
+
+        return 0;
+    }
+
     Scalar getMicrobialConcentration(unsigned elemIdx) const
     {
         if (cMicrobes_.size() > elemIdx)
@@ -204,6 +220,7 @@ protected:
                                    bool enablePolymer,
                                    bool enableFoam,
                                    bool enableBrine,
+                                   bool enableSaltPrecipitation,
                                    bool enableExtbo,
                                    bool enableMICP);
 
@@ -360,6 +377,7 @@ protected:
     bool enablePolymer_;
     bool enableFoam_;
     bool enableBrine_;
+    bool enableSaltPrecipitation_;
     bool enableExtbo_;
     bool enableMICP_;
 
@@ -393,6 +411,8 @@ protected:
     ScalarBuffer cPolymer_;
     ScalarBuffer cFoam_;
     ScalarBuffer cSalt_;
+    ScalarBuffer pSalt_;
+    ScalarBuffer permFact_;
     ScalarBuffer extboX_;
     ScalarBuffer extboY_;
     ScalarBuffer extboZ_;

ebos/ecloutputblackoilmodule.hh

@@ -135,6 +135,7 @@ public:
                    getPropValue<TypeTag, Properties::EnablePolymer>(),
                    getPropValue<TypeTag, Properties::EnableFoam>(),
                    getPropValue<TypeTag, Properties::EnableBrine>(),
+                   getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
                    getPropValue<TypeTag, Properties::EnableExtbo>(),
                    getPropValue<TypeTag, Properties::EnableMICP>())
          , simulator_(simulator)
@@ -332,6 +333,14 @@ public:
                 this->cSalt_[globalDofIdx] = fs.saltConcentration().value();
             }
 
+            if (!this->pSalt_.empty()) {
+                this->pSalt_[globalDofIdx] = intQuants.saltSaturation().value();
+            }
+
+            if (!this->permFact_.empty()) {
+                this->permFact_[globalDofIdx] = intQuants.permFactor().value();
+            }
+            
             if (!this->extboX_.empty()) {
                 this->extboX_[globalDofIdx] = intQuants.xVolume().value();
             }

ebos/eclproblem.hh

@@ -619,6 +619,7 @@ class EclProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
     enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
+    enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
     enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
@@ -1830,8 +1831,14 @@ public:
         if constexpr (enablePolymerMolarWeight)
             values[Indices::polymerMoleWeightIdx]= this->polymerMoleWeight_[globalDofIdx];
 
-        if constexpr (enableBrine)
+        if constexpr (enableBrine) {
-            values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
+            if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::Sp) {
+                values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation();
+            }
+            else {
+                values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration();
+            }
+        }
 
         if constexpr (enableMICP){
             values[Indices::microbialConcentrationIdx]= this->microbialConcentration_[globalDofIdx];

ebos/equil/initstateequil.hh

@@ -47,6 +47,7 @@
 #include <opm/input/eclipse/EclipseState/Tables/PbvdTable.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/PdvdTable.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/SaltvdTable.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/SaltpvdTable.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <fmt/format.h>
 
@@ -1580,6 +1581,7 @@ public:
                          const bool applySwatInit = true)
         : temperature_(gridView.size(/*codim=*/0)),
           saltConcentration_(gridView.size(/*codim=*/0)),
+          saltSaturation_(gridView.size(/*codim=*/0)),
           pp_(FluidSystem::numPhases,
           std::vector<double>(gridView.size(/*codim=*/0))),
           sat_(FluidSystem::numPhases,
@@ -1711,6 +1713,9 @@ public:
         // EXTRACT the initial salt concentration
         updateInitialSaltConcentration_(eclipseState, eqlmap);
 
+        // EXTRACT the initial salt saturation
+        updateInitialSaltSaturation_(eclipseState, eqlmap);
+
         // Compute pressures, saturations, rs and rv factors.
         const auto& comm = gridView.comm();
         calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
@@ -1727,6 +1732,7 @@ public:
 
     const Vec& temperature() const { return temperature_; }
     const Vec& saltConcentration() const { return saltConcentration_; }
+    const Vec& saltSaturation() const { return saltSaturation_; }
     const PVec& press() const { return pp_; }
     const PVec& saturation() const { return sat_; }
     const Vec& rs() const { return rs_; }
@@ -1767,14 +1773,36 @@ private:
             }
         }
     }
+    template <class RMap>
+    void updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg)
+    {
+        const int numEquilReg = rsFunc_.size();
+        saltpVdTable_.resize(numEquilReg);
+        const auto& tables = eclState.getTableManager();
+        const TableContainer& saltpvdTables = tables.getSaltpvdTables();
+
+        for (size_t i = 0; i < saltpvdTables.size(); ++i) {
+            const SaltpvdTable& saltpvdTable = saltpvdTables.getTable<SaltpvdTable>(i);
+            saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn());
+
+            const auto& cells = reg.cells(i);
+            for (const auto& cell : cells) {
+                const double depth = cellCenterDepth_[cell];
+                this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true);
+            }
+        }
+    }
+
 
     std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
     std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
     using TabulatedFunction = Tabulated1DFunction<double>;
     std::vector<TabulatedFunction> saltVdTable_;
+    std::vector<TabulatedFunction> saltpVdTable_;
     std::vector<int> regionPvtIdx_;
     Vec temperature_;
     Vec saltConcentration_;
+    Vec saltSaturation_;
     PVec pp_;
     PVec sat_;
     Vec rs_;

flow/flow_brine_saltprecipitation.cpp

@@ -0,0 +1,24 @@
+/*
+  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"
+#include <flow/flow_ebos_brine_saltprecipitation.hpp>
+
+
+int main(int argc, char** argv)
+{
+    return Opm::flowEbosBrineSaltPrecipitationMainStandalone(argc, argv);
+}

flow/flow_ebos_brine_saltprecipitation.cpp

@@ -0,0 +1,80 @@
+/*
+  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"
+
+#include <flow/flow_ebos_brine_saltprecipitation.hpp>
+
+#include <opm/material/common/ResetLocale.hpp>
+#include <opm/grid/CpGrid.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/simulators/flow/Main.hpp>
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowBrineSaltPrecipitationProblem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+template<class TypeTag>
+struct EnableBrine<TypeTag, TTag::EclFlowBrineSaltPrecipitationProblem> {
+    static constexpr bool value = true;
+};
+
+template<class TypeTag>
+struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowBrineSaltPrecipitationProblem> {
+    static constexpr bool value = true;
+};
+}}
+
+namespace Opm {
+void flowEbosBrineSaltPrecipitationSetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                          std::shared_ptr<EclipseState> eclState,
+                          std::shared_ptr<Schedule> schedule,
+                          std::shared_ptr<SummaryConfig> summaryConfig)
+{
+    using TypeTag = Properties::TTag::EclFlowBrineSaltPrecipitationProblem;
+    using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
+
+    Vanguard::setExternalSetupTime(setupTime);
+    Vanguard::setExternalDeck(std::move(deck));
+    Vanguard::setExternalEclState(std::move(eclState));
+    Vanguard::setExternalSchedule(std::move(schedule));
+    Vanguard::setExternalSummaryConfig(std::move(summaryConfig));
+}
+
+
+// ----------------- Main program -----------------
+int flowEbosBrineSaltPrecipitationMain(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMainEbos<Properties::TTag::EclFlowBrineSaltPrecipitationProblem>
+        mainfunc {argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+int flowEbosBrineSaltPrecipitationMainStandalone(int argc, char** argv)
+{
+    using TypeTag = Properties::TTag::EclFlowBrineSaltPrecipitationProblem;
+    auto mainObject = Opm::Main(argc, argv);
+    return mainObject.runStatic<TypeTag>();
+}
+
+}

flow/flow_ebos_brine_saltprecipitation.hpp

@@ -0,0 +1,42 @@
+/*
+  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/>.
+*/
+#ifndef FLOW_EBOS_BRINE_SALTPRECIPITATION_HPP
+#define FLOW_EBOS_BRINE_SALTPRECIPITATION_HPP
+
+#include <memory>
+
+namespace Opm {
+
+class Deck;
+class EclipseState;
+class Schedule;
+class SummaryConfig;
+
+void flowEbosBrineSaltPrecipitationSetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                          std::shared_ptr<EclipseState> eclState,
+                          std::shared_ptr<Schedule> schedule,
+                          std::shared_ptr<SummaryConfig> summaryConfig);
+
+//! \brief Main function used in flow binary.
+int flowEbosBrineSaltPrecipitationMain(int argc, char** argv, bool outputCout, bool outputFiles);
+
+//! \brief Main function used in flow_brine binary.
+int flowEbosBrineSaltPrecipitationMainStandalone(int argc, char** argv);
+
+}
+
+#endif // FLOW_EBOS_BRINE_HPP

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -127,6 +127,10 @@ struct EnableBrine<TypeTag, TTag::EclFlowProblem> {
     static constexpr bool value = false;
 };
 template<class TypeTag>
+struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
 struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
     static constexpr bool value = false;
 };

opm/simulators/flow/Main.hpp

@@ -33,6 +33,7 @@
 #include <flow/flow_ebos_extbo.hpp>
 #include <flow/flow_ebos_foam.hpp>
 #include <flow/flow_ebos_brine.hpp>
+#include <flow/flow_ebos_brine_saltprecipitation.hpp>
 #include <flow/flow_ebos_oilwater_brine.hpp>
 #include <flow/flow_ebos_gaswater_brine.hpp>
 #include <flow/flow_ebos_energy.hpp>
@@ -668,6 +669,11 @@ private:
                 return flowEbosGasWaterBrineMain(argc_, argv_, outputCout_, outputFiles_);
             }
         }
+        else if (eclipseState_->getSimulationConfig().hasPRECSALT()) {
+            flowEbosBrineSaltPrecipitationSetDeck(
+                setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
+            return flowEbosBrineSaltPrecipitationMain(argc_, argv_, outputCout_, outputFiles_);
+        }
         else {
             flowEbosBrineSetDeck(
                 setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -472,7 +472,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"PEKTABX", {false, std::nullopt}},
         {"PENUM", {false, std::nullopt}},
         {"PERMAVE", {false, std::nullopt}},
-        {"PERMFACT", {false, std::nullopt}},
         {"PERMXY", {false, std::nullopt}},
         {"PERMYZ", {false, std::nullopt}},
         {"PERMZX", {false, std::nullopt}},
@@ -488,7 +487,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"PLYCAMAX", {false, std::nullopt}},
         {"PLYDHFLF", {false, std::nullopt}},
         {"PPCWMAX", {false, std::nullopt}},
-        {"PRECSALT", {false, std::nullopt}},
         {"PRORDER", {false, std::nullopt}},
         {"PRVD", {false, std::nullopt}},
         {"PVTGW", {false, std::nullopt}},
@@ -553,8 +551,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"RSSPEC", {false, std::nullopt}},
         {"RTEMPA", {false, std::nullopt}},
         {"RWGSALT", {false, std::nullopt}},
-        {"SALTPVD", {false, std::nullopt}},
-        {"SALTSOL", {false, std::nullopt}},
         {"SAMG", {false, std::nullopt}},
         {"SAVE", {false, std::nullopt}},
         {"SKIP", {false, std::nullopt}},