blackoilbrinemodules.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  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 2 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/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_BLACK_OIL_BRINE_MODULE_HH
#define EWOMS_BLACK_OIL_BRINE_MODULE_HH
 
#include "blackoilproperties.hh"
 
#include <opm/models/blackoil/blackoilbrineparams.hh>
 
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PvtwsaltTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PcfactTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PermfactTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SaltSolubilityTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#endif
 
#include <dune/common/fvector.hh>
 
#include <string>
#include <math.h>
 
namespace Opm {
template <class TypeTag, bool enableBrineV = getPropValue<TypeTag, Properties::EnableBrine>()>
class BlackOilBrineModule
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using Model = GetPropType<TypeTag, Properties::Model>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
 
    using Toolbox = MathToolbox<Evaluation>;
 
    using TabulatedFunction = typename BlackOilBrineParams<Scalar>::TabulatedFunction;
 
    static constexpr unsigned saltConcentrationIdx = Indices::saltConcentrationIdx;
    static constexpr unsigned contiBrineEqIdx = Indices::contiBrineEqIdx;
    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
    static constexpr bool gasEnabled = Indices::gasEnabled;
    static constexpr bool oilEnabled = Indices::oilEnabled;
    static constexpr unsigned enableBrine = enableBrineV;
    static constexpr unsigned enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
 
    static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
    static constexpr unsigned numPhases = FluidSystem::numPhases;
 
public:
 
#if HAVE_ECL_INPUT
    static void initFromState(const EclipseState& eclState)
    {
        // some sanity checks: if brine are enabled, the BRINE keyword must be
        // present, if brine are disabled the keyword must not be present.
        if (enableBrine && !eclState.runspec().phases().active(Phase::BRINE)) {
            throw std::runtime_error("Non-trivial brine treatment requested at compile time, but "
                                     "the deck does not contain the BRINE keyword");
        }
        else if (!enableBrine && eclState.runspec().phases().active(Phase::BRINE)) {
            throw std::runtime_error("Brine treatment disabled at compile time, but the deck "
                                     "contains the BRINE keyword");
        }
 
        if (!eclState.runspec().phases().active(Phase::BRINE))
            return; // brine treatment is supposed to be disabled
 
        const auto& tableManager = eclState.getTableManager();
 
        unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables();
        params_.referencePressure_.resize(numPvtRegions);
 
        const auto& pvtwsaltTables = tableManager.getPvtwSaltTables();
 
        // initialize the objects which deal with the BDENSITY keyword
        const auto& bdensityTables = tableManager.getBrineDensityTables();
        if (!bdensityTables.empty()) {
            params_.bdensityTable_.resize(numPvtRegions);
            assert(numPvtRegions == bdensityTables.size());
            for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
                const auto& bdensityTable = bdensityTables[pvtRegionIdx];
                const auto& pvtwsaltTable = pvtwsaltTables[pvtRegionIdx];
                const auto& c = pvtwsaltTable.getSaltConcentrationColumn();
                params_.bdensityTable_[pvtRegionIdx].setXYContainers(c, bdensityTable);
            }
        }
 
        if constexpr (enableSaltPrecipitation) {
            const TableContainer& permfactTables = tableManager.getPermfactTables();
            params_.permfactTable_.resize(numPvtRegions);
            for (size_t i = 0; i < permfactTables.size(); ++i) {
                const PermfactTable& permfactTable = permfactTables.getTable<PermfactTable>(i);
                params_.permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn());
            }
 
            const TableContainer& saltsolTables = tableManager.getSaltsolTables();
            if (!saltsolTables.empty()) {
                params_.saltsolTable_.resize(numPvtRegions);
                params_.saltdenTable_.resize(numPvtRegions);
                assert(numPvtRegions == saltsolTables.size());
                for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
                    const SaltsolTable& saltsolTable = saltsolTables.getTable<SaltsolTable>(pvtRegionIdx );
                    params_.saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front();
                    params_.saltdenTable_[pvtRegionIdx] = saltsolTable.getSaltdenColumn().front();
                }
            }
 
            const TableContainer& pcfactTables = tableManager.getPcfactTables();
            if (!pcfactTables.empty()) {
                unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
                params_.pcfactTable_.resize(numSatRegions);
                for (size_t i = 0; i < pcfactTables.size(); ++i) {
                    const PcfactTable& pcfactTable = pcfactTables.getTable<PcfactTable>(i);
                    params_.pcfactTable_[i].setXYContainers(pcfactTable.getPorosityChangeColumn(), pcfactTable.getPcMultiplierColumn());
                }
            }
        }
    }
#endif
 
    static void registerParameters()
    {
    }
 
    static bool primaryVarApplies(unsigned pvIdx)
    {
        if constexpr (enableBrine)
            return pvIdx == saltConcentrationIdx;
        else
            return false;
    }
 
    template <class FluidState>
    static void assignPrimaryVars(PrimaryVariables& priVars,
                                  const FluidState& fluidState)
    {
        if constexpr (enableBrine)
            priVars[saltConcentrationIdx] = fluidState.saltConcentration();
    }
 
    static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
    {
        assert(primaryVarApplies(pvIdx));
 
        return "saltConcentration";
    }
 
    static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
    {
        assert(primaryVarApplies(pvIdx));
 
        // TODO: it may be beneficial to chose this differently.
        return static_cast<Scalar>(1.0);
    }
 
    static bool eqApplies(unsigned eqIdx)
    {
        if constexpr (enableBrine)
            return eqIdx == contiBrineEqIdx;
        else
            return false;
    }
 
    static std::string eqName([[maybe_unused]] unsigned eqIdx)
    {
        assert(eqApplies(eqIdx));
 
        return "conti^brine";
    }
 
    static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
    {
        assert(eqApplies(eqIdx));
 
        // TODO: it may be beneficial to chose this differently.
        return static_cast<Scalar>(1.0);
    }
 
    // must be called after water storage is computed
    template <class LhsEval>
    static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                           const IntensiveQuantities& intQuants)
    {
        if constexpr (enableBrine) {
            const auto& fs = intQuants.fluidState();
 
            LhsEval surfaceVolumeWater =
                    Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
                    * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
                    * Toolbox::template decay<LhsEval>(intQuants.porosity());
 
            // avoid singular matrix if no water is present.
            surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
 
            // Brine in water phase
            const LhsEval massBrine = surfaceVolumeWater
                    * Toolbox::template decay<LhsEval>(fs.saltConcentration());
 
            if (enableSaltPrecipitation){
                double saltDensity = intQuants.saltDensity(); // Solid salt density kg/m3
                const LhsEval solidSalt =
                              Toolbox::template decay<LhsEval>(intQuants.porosity())
                              / (1.0 - Toolbox::template decay<LhsEval>(intQuants.saltSaturation()) + 1.e-8)
                              * saltDensity
                              * Toolbox::template decay<LhsEval>(intQuants.saltSaturation());
 
                storage[contiBrineEqIdx] += massBrine + solidSalt;
            }
            else { storage[contiBrineEqIdx] += massBrine;}
        }
    }
 
    static void computeFlux([[maybe_unused]] RateVector& flux,
                            [[maybe_unused]] const ElementContext& elemCtx,
                            [[maybe_unused]] unsigned scvfIdx,
                            [[maybe_unused]] unsigned timeIdx)
 
    {
        if constexpr (enableBrine) {
            const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
 
            const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
            const unsigned inIdx = extQuants.interiorIndex();
            const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
 
            if (upIdx == inIdx) {
                flux[contiBrineEqIdx] =
                        extQuants.volumeFlux(waterPhaseIdx)
                        *up.fluidState().invB(waterPhaseIdx)
                        *up.fluidState().saltConcentration();
            }
            else {
                flux[contiBrineEqIdx] =
                        extQuants.volumeFlux(waterPhaseIdx)
                        *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
                        *decay<Scalar>(up.fluidState().saltConcentration());
            }
        }
    }
 
    static Scalar computeUpdateError(const PrimaryVariables&,
                                     const EqVector&)
    {
        // do not consider consider the change of Brine primary variables for
        // convergence
        // TODO: maybe this should be changed
        return static_cast<Scalar>(0.0);
    }
 
    template <class DofEntity>
    static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
    {
        if constexpr (enableBrine) {
            unsigned dofIdx = model.dofMapper().index(dof);
            const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
            outstream << priVars[saltConcentrationIdx];
        }
    }
 
    template <class DofEntity>
    static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
    {
        if constexpr (enableBrine) {
            unsigned dofIdx = model.dofMapper().index(dof);
            PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
            PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
 
            instream >> priVars0[saltConcentrationIdx];
 
            // set the primary variables for the beginning of the current time step.
            priVars1[saltConcentrationIdx] = priVars0[saltConcentrationIdx];
        }
    }
 
    static const Scalar& referencePressure(const ElementContext& elemCtx,
                                           unsigned scvIdx,
                                           unsigned timeIdx)
    {
        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
        return params_.referencePressure_[pvtnumRegionIdx];
    }
 
 
    static const TabulatedFunction& bdensityTable(const ElementContext& elemCtx,
                                                  unsigned scvIdx,
                                                  unsigned timeIdx)
    {
        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
        return params_.bdensityTable_[pvtnumRegionIdx];
    }
 
    static const TabulatedFunction& pcfactTable(unsigned satnumRegionIdx)
    {
        return params_.pcfactTable_[satnumRegionIdx];
    }
 
    static const TabulatedFunction& permfactTable(const ElementContext& elemCtx,
                                                  unsigned scvIdx,
                                                  unsigned timeIdx)
    {
        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
        return params_.permfactTable_[pvtnumRegionIdx];
    }
 
    static const TabulatedFunction& permfactTable(unsigned pvtnumRegionIdx)
    {
        return params_.permfactTable_[pvtnumRegionIdx];
    }
 
    static const Scalar saltsolTable(const ElementContext& elemCtx,
                                                  unsigned scvIdx,
                                                  unsigned timeIdx)
    {
        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
        return params_.saltsolTable_[pvtnumRegionIdx];
    }
 
    static const Scalar saltdenTable(const ElementContext& elemCtx,
                                                  unsigned scvIdx,
                                                  unsigned timeIdx)
    {
        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
        return params_.saltdenTable_[pvtnumRegionIdx];
    }
 
    static bool hasBDensityTables()
    {
        return !params_.bdensityTable_.empty();
    }
 
    static bool hasSaltsolTables()
    {
        return !params_.saltsolTable_.empty();
    }
 
    static bool hasPcfactTables()
    {
        if constexpr (enableSaltPrecipitation)
            return !params_.pcfactTable_.empty();
        else
            return false;
    }
 
    static Scalar saltSol(unsigned regionIdx) {
        return params_.saltsolTable_[regionIdx];
    }
 
private:
    static BlackOilBrineParams<Scalar> params_;
};
 
 
template <class TypeTag, bool enableBrineV>
BlackOilBrineParams<typename BlackOilBrineModule<TypeTag, enableBrineV>::Scalar>
BlackOilBrineModule<TypeTag, enableBrineV>::params_;
 
template <class TypeTag, bool enableBrineV = getPropValue<TypeTag, Properties::EnableBrine>()>
class BlackOilBrineIntensiveQuantities
{
    using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
 
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
 
    using BrineModule = BlackOilBrineModule<TypeTag>;
 
    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
    static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
    static constexpr unsigned enableBrine = enableBrineV;
    static constexpr unsigned enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
    static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
 
public:
 
    void updateSaltConcentration_(const ElementContext& elemCtx,
                                  unsigned dofIdx,
                                  unsigned timeIdx)
    {
        const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
 
        auto& fs = asImp_().fluidState_;
 
        if constexpr (enableSaltPrecipitation) {
            const auto& saltsolTable = BrineModule::saltsolTable(elemCtx, dofIdx, timeIdx);
            saltSolubility_ = saltsolTable;
 
            const auto& saltdenTable = BrineModule::saltdenTable(elemCtx, dofIdx, timeIdx);
            saltDensity_ = saltdenTable;
 
            if (priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
                saltSaturation_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
                fs.setSaltConcentration(saltSolubility_);
            }
            else {
                saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
                fs.setSaltConcentration(saltConcentration_);
                saltSaturation_ = 0.0;
            }
            fs.setSaltSaturation(saltSaturation_);
        }
        else {
            saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
            fs.setSaltConcentration(saltConcentration_);
        }
    }
 
    void saltPropertiesUpdate_([[maybe_unused]] const ElementContext& elemCtx,
                               [[maybe_unused]] unsigned dofIdx,
                               [[maybe_unused]] unsigned timeIdx)
    {
        if constexpr (enableSaltPrecipitation) {
            const Evaluation porosityFactor  = min(1.0 - saltSaturation(), 1.0); //phi/phi_0
 
            const auto& permfactTable = BrineModule::permfactTable(elemCtx, dofIdx, timeIdx);
 
            permFactor_ = permfactTable.eval(porosityFactor);
            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                if (!FluidSystem::phaseIsActive(phaseIdx))
                    continue;
 
                asImp_().mobility_[phaseIdx] *= permFactor_;
            }
        }
    }
 
    const Evaluation& saltConcentration() const
    { return saltConcentration_; }
 
    const Evaluation& brineRefDensity() const
    { return refDensity_; }
 
    const Evaluation& saltSaturation() const
    { return saltSaturation_; }
 
    Scalar saltSolubility() const
    { return saltSolubility_; }
 
    Scalar saltDensity() const
    { return saltDensity_; }
 
    const Evaluation& permFactor() const
    { return permFactor_; }
 
protected:
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }
 
    Evaluation saltConcentration_;
    Evaluation refDensity_;
    Evaluation saltSaturation_;
    Evaluation permFactor_;
    Scalar saltSolubility_;
    Scalar saltDensity_;
 
};
 
template <class TypeTag>
class BlackOilBrineIntensiveQuantities<TypeTag, false>
{
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
public:
    void updateSaltConcentration_(const ElementContext&,
                                  unsigned,
                                  unsigned)
    { }
 
    void saltPropertiesUpdate_(const ElementContext&,
                                  unsigned,
                                  unsigned)
    { }
 
    const Evaluation& saltConcentration() const
    { throw std::runtime_error("saltConcentration() called but brine are disabled"); }
 
    const Evaluation& brineRefDensity() const
    { throw std::runtime_error("brineRefDensity() called but brine are disabled"); }
 
    const Evaluation& saltSaturation() const
    { throw std::logic_error("saltSaturation() called but salt precipitation is disabled"); }
 
    const Scalar saltSolubility() const
    { throw std::logic_error("saltSolubility() called but salt precipitation is disabled"); }
 
    const Scalar saltDensity() const
    { throw std::logic_error("saltDensity() called but salt precipitation is disabled"); }
 
    const Evaluation& permFactor() const
    { throw std::logic_error("permFactor() called but salt precipitation is disabled"); }
 
};
 
} // namespace Ewoms
 
#endif