diff --git a/opm/models/blackoil/blackoildispersionmodule.hh b/opm/models/blackoil/blackoildispersionmodule.hh new file mode 100644 index 000000000..640102048 --- /dev/null +++ b/opm/models/blackoil/blackoildispersionmodule.hh @@ -0,0 +1,520 @@ +// -*- 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 . + + 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. +*/ +/*! + * \file + * + * \brief Classes required for mechanical dispersion. + */ +#ifndef EWOMS_DISPERSION_MODULE_HH +#define EWOMS_DISPERSION_MODULE_HH + +#include + +#include + +#include + +#include + +namespace Opm { + +/*! + * \ingroup Dispersion + * \class Opm::BlackOilDispersionModule + * \brief Provides the auxiliary methods required for consideration of the + * dispersion equation. + */ +template +class BlackOilDispersionModule; + +template +class BlackOilDispersionExtensiveQuantities; + +/*! + * \copydoc Opm::BlackOilDispersionModule + */ +template +class BlackOilDispersionModule +{ + using Scalar = GetPropType; + using RateVector = GetPropType; + using FluidSystem = GetPropType; + using Evaluation = GetPropType; + + enum { numPhases = FluidSystem::numPhases }; + +public: + using ExtensiveQuantities = BlackOilDispersionExtensiveQuantities; + /*! + * \brief Adds the dispersive flux to the flux vector over a flux + * integration point. + */ + template + static void addDispersiveFlux(RateVector&, + const Context&, + unsigned, + unsigned) + {} + + template + static void addDispersiveFlux(RateVector&, + const FluidState&, + const FluidState&, + const Evaluation&, + const Scalar&) + {} +}; + +/*! + * \copydoc Opm::BlackOilDispersionModule + */ +template +class BlackOilDispersionModule +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using PrimaryVariables = GetPropType; + using IntensiveQuantities = GetPropType; + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + using Model = GetPropType; + using Simulator = GetPropType; + using EqVector = GetPropType; + using RateVector = GetPropType; + using Indices = GetPropType; + + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { enableDispersion = getPropValue() }; + + using Toolbox = MathToolbox; + +public: + using ExtensiveQuantities = BlackOilDispersionExtensiveQuantities; + /*! + * \brief Adds the mass flux due to dispersion to the flux vector over the + * flux integration point. + */ + template + static void addDispersiveFlux(RateVector& flux, const Context& context, + unsigned spaceIdx, unsigned timeIdx) + { + // Only work if dispersion is enabled by DISPERC in the deck + if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) { + return; + } + const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); + const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState(); + const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState(); + const auto& dispersivity = extQuants.dispersivity(); + const auto& normVelocityAvg = extQuants.normVelocityAvg(); + addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg); + } + + /*! + * \brief Adds the mass flux due to dispersion to the flux vector over the + * integration point. Following the notation in blackoilmodel.hh, + * the dispersive flux for component \f$\kappa\f$ in phase \f$\alpha\f$ + * is given by: \f$-b_\alpha E||\mathrm{v}_\alpha||\mathbf{grad}X_\alpha^\kappa\f$, + * where \f$b_\alpha\f$ is the shrinkage/expansion factor [-], E is the isotropic + * dispersivity coefficient [L], \f$\mathrm{v}_\alpha\f$ is the filter velocity + * [L/T], and \f$X_\alpha^\kappa\f$ the component mass fraction [-]. Each component mass + * fraction can be computed using \f$R_s,\;R_v,\;R_{sw},\;R_{vw}\f$. For example, + * \f$X_w^G=\frac{R_{sw}}{R_{sw}+\rho_w/\rho_g}\f$, where \f$\rho_w\f$ and \f$\rho_g\f$ + * are the reference densities. + * Following the implementation of the diffusive flux (blackoildiffusionmodule.hh) and considering + * the case for the water phase and gas component as an example, for cells i and j, the discrete version + * of the dispersive flux at the face's integration point is given by + * \f$-b_{w,ij}v_{w,ij}(\frac{1}{R_{sw,ij}+\rho_w/\rho_g})D_{ij}(R_{sw,i}-R_{sw,j})\f$ + * where \f$b_{w,ij}\f$, \f$v_{w,ij}\f$, and \f$R_{sw,ij}\f$ are computed using the arithmetic mean, and + * the ratio \f$\frac{1}{R_{sw,ij}+\rho_w/\rho_g}\f$ is denoted as conversion factor. The dispersivity + * \f$D_{ij}\f$ is computed in ecltransmissibility_impl.hh, using the dispersion coefficients \f$E_i\f$ + * and \f$E_j\f$. + */ + template + static void addDispersiveFlux(RateVector& flux, + const FluidState& fluidStateI, + const FluidState& fluidStateJ, + const Evaluation& dispersivity, + const Scalar& normVelocityAvg) + { + unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + // no dispersion in water for blackoil models unless water can contain dissolved gas + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + + // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil + if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) { + continue; + } + + // arithmetic mean of the phase's b factor + Evaluation bAvg = fluidStateI.invB(phaseIdx); + bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx)); + bAvg /= 2; + + Evaluation convFactor = 1.0; + Evaluation diffR = 0.0; + if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { + Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2; + convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg); + diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs()); + } + if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2; + convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex)); + diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv()); + } + if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { + Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2; + convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg); + diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw()); + } + if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2; + convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex)); + diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw()); + } + + // mass flux of solvent component + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + flux[conti0EqIdx + activeSolventCompIdx] += + - bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; + + // mass flux of solute component + unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); + unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); + flux[conti0EqIdx + activeSoluteCompIdx] += + bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; + } + } + +private: + static Scalar toMassFractionGasOil (unsigned regionIdx) { + Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx); + Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx); + return rhoO / rhoG; + } + static Scalar toMassFractionGasWater (unsigned regionIdx) { + Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx); + Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx); + return rhoW / rhoG; + } +}; + +/*! + * \ingroup Dispersion + * \class Opm::BlackOilDispersionIntensiveQuantities + * + * \brief Provides the volumetric quantities required for the + * calculation of dispersive fluxes. + */ +template +class BlackOilDispersionIntensiveQuantities; + +/*! + * \copydoc Opm::DispersionIntensiveQuantities + */ +template +class BlackOilDispersionIntensiveQuantities +{ + using Scalar = GetPropType; + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + +public: + /*! + * \brief Returns the max. norm of the filter velocity of the cell. + */ + Scalar normVelocityCell(unsigned, unsigned) const + { + throw std::logic_error("Method normVelocityCell() " + "does not make sense if dispersion is disabled"); + } + +protected: + /*! + * \brief Update the quantities required to calculate dispersive + * fluxes. + */ + template + void update_(ElementContext&, + unsigned, + unsigned) + { } +}; + +/*! + * \copydoc Opm::DispersionIntensiveQuantities + */ +template +class BlackOilDispersionIntensiveQuantities +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + using IntensiveQuantities = GetPropType; + using Indices = GetPropType; + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + enum { gasCompIdx = FluidSystem::gasCompIdx }; + enum { oilCompIdx = FluidSystem::oilCompIdx }; + enum { waterCompIdx = FluidSystem::waterCompIdx }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { enableDispersion = getPropValue() }; + +public: + /*! + * \brief Returns the max. norm of the filter velocity of the cell. + */ + Scalar normVelocityCell(unsigned phaseIdx) const + { + + return normVelocityCell_[phaseIdx]; + } + +protected: + /*! + * \brief Update the quantities required to calculate dispersive + * mass fluxes. This considers the linear disperison model + * described in the SPE CSP11 benchmark document (eq. 2.3) + * https://github.com/Simulation-Benchmarks/11thSPE-CSP/ + * blob/main/description/spe_csp11_description.pdf + * The maximum norm is used to compute the cell + * filter velocity value of the corresponding phase. + */ + template + void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + { + // Only work if dispersion is enabled by DISPERC in the deck + if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) { + return; + } + const auto& problem = elemCtx.simulator().problem(); + if (problem.model().linearizer().getVelocityInfo().empty()) { + return; + } + const std::array phaseIdxs = { gasPhaseIdx, oilPhaseIdx, waterPhaseIdx }; + const std::array compIdxs = { gasCompIdx, oilCompIdx, waterCompIdx }; + const auto& velocityInf = problem.model().linearizer().getVelocityInfo(); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); + auto velocityInfos = velocityInf[globalDofIdx]; + for (unsigned i = 0; i < phaseIdxs.size(); ++i) { + normVelocityCell_[i] = 0; + } + for (auto& velocityInfo : velocityInfos) { + for (unsigned i = 0; i < phaseIdxs.size(); ++i) { + if (FluidSystem::phaseIsActive(phaseIdxs[i])) { + normVelocityCell_[phaseIdxs[i]] = max( normVelocityCell_[phaseIdxs[i]], + std::abs( velocityInfo.velocity[conti0EqIdx + + Indices::canonicalToActiveComponentIndex(compIdxs[i])] )); + } + } + } + } + +private: + Scalar normVelocityCell_[numPhases]; +}; + +/*! + * \ingroup Dispersion + * \class Opm::BlackOilDispersionExtensiveQuantities + * + * \brief Provides the quantities required to calculate dispersive mass fluxes. + */ +template +class BlackOilDispersionExtensiveQuantities; + +/*! + * \copydoc Opm::DispersionExtensiveQuantities + */ +template +class BlackOilDispersionExtensiveQuantities +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + using IntensiveQuantities = GetPropType; + + enum { numPhases = FluidSystem::numPhases }; + +protected: + /*! + * \brief Update the quantities required to calculate + * the dispersive fluxes. + */ + void update_(const ElementContext&, + unsigned, + unsigned) + {} + + template + void updateBoundary_(const Context&, + unsigned, + unsigned, + const FluidState&) + {} + +public: + using ScalarArray = Scalar[numPhases]; + + static void update(ScalarArray&, + const IntensiveQuantities&, + const IntensiveQuantities&) + {} + + /*! + * \brief The dispersivity the face. + * + */ + const Scalar& dispersivity() const + { + throw std::logic_error("The method dispersivity() does not " + "make sense if dispersion is disabled."); + } + + /*! + * \brief The effective filter velocity coefficient in a + * fluid phase at the face's integration point + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Scalar& normVelocityAvg(unsigned) const + { + throw std::logic_error("The method normVelocityAvg() " + "does not make sense if dispersion is disabled."); + } + +}; + +/*! + * \copydoc Opm::BlackOilDispersionExtensiveQuantities + */ +template +class BlackOilDispersionExtensiveQuantities +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using ElementContext = GetPropType; + using GridView = GetPropType; + using FluidSystem = GetPropType; + using Toolbox = MathToolbox; + using IntensiveQuantities = GetPropType; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + + using DimVector = Dune::FieldVector; + using DimEvalVector = Dune::FieldVector; +public: + using ScalarArray = Scalar[numPhases]; + static void update(ScalarArray& normVelocityAvg, + const IntensiveQuantities& intQuantsInside, + const IntensiveQuantities& intQuantsOutside) + { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + // no dispersion in water for blackoil models unless water can contain dissolved gas + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil + if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) { + continue; + } + // use the arithmetic average for the effective + // velocity coefficients at the face's integration point. + normVelocityAvg[phaseIdx] = 0.5 * + ( intQuantsInside.normVelocityCell(phaseIdx) + + intQuantsOutside.normVelocityCell(phaseIdx) ); + Valgrind::CheckDefined(normVelocityAvg[phaseIdx]); + } + } +protected: + template + void updateBoundary_(const Context&, + unsigned, + unsigned, + const FluidState&) + { + throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil"); + } + +public: + /*! + * \brief The dispersivity of the face. + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Scalar& dispersivity() const + { return dispersivity_; } + + /*! + * \brief The effective velocity coefficient in a + * fluid phase at the face's integration point + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Scalar& normVelocityAvg(unsigned phaseIdx) const + { return normVelocityAvg_[phaseIdx]; } + + const auto& normVelocityAvg() const{ + return normVelocityAvg_; + } + +private: + Scalar dispersivity_; + ScalarArray normVelocityAvg_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index db00dc53f..8a4e19878 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -36,6 +36,7 @@ #include "blackoilbrinemodules.hh" #include "blackoilenergymodules.hh" #include "blackoildiffusionmodule.hh" +#include "blackoildispersionmodule.hh" #include "blackoilmicpmodules.hh" #include @@ -70,6 +71,7 @@ class BlackOilIntensiveQuantities : public GetPropType , public GetPropType::FluxIntensiveQuantities , public BlackOilDiffusionIntensiveQuantities() > + , public BlackOilDispersionIntensiveQuantities() > , public BlackOilSolventIntensiveQuantities , public BlackOilExtboIntensiveQuantities , public BlackOilPolymerIntensiveQuantities @@ -103,6 +105,7 @@ class BlackOilIntensiveQuantities enum { enableTemperature = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; + enum { enableDispersion = getPropValue() }; enum { enableMICP = getPropValue() }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; @@ -124,6 +127,7 @@ class BlackOilIntensiveQuantities using DimMatrix = Dune::FieldMatrix; using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities; using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities; + using DispersionIntensiveQuantities = BlackOilDispersionIntensiveQuantities; using DirectionalMobilityPtr = Opm::Utility::CopyablePtr>; @@ -473,6 +477,9 @@ public: // update the diffusion specific quantities of the intensive quantities DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx); + // update the dispersion specific quantities of the intensive quantities + DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx); + #ifndef NDEBUG // some safety checks in debug mode for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index abf267e83..e3a058af7 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -36,6 +36,7 @@ #include "blackoilfoammodules.hh" #include "blackoilbrinemodules.hh" #include "blackoildiffusionmodule.hh" +#include "blackoildispersionmodule.hh" #include "blackoilmicpmodules.hh" #include #include @@ -92,6 +93,7 @@ class BlackOilLocalResidualTPFA : public GetPropType(); static constexpr bool enableBrine = getPropValue(); static constexpr bool enableDiffusion = getPropValue(); + static constexpr bool enableDispersion = getPropValue(); static constexpr bool enableMICP = getPropValue(); using SolventModule = BlackOilSolventModule; @@ -101,6 +103,7 @@ class BlackOilLocalResidualTPFA : public GetPropType; using BrineModule = BlackOilBrineModule; using DiffusionModule = BlackOilDiffusionModule; + using DispersionModule = BlackOilDispersionModule; using MICPModule = BlackOilMICPModule; using Toolbox = MathToolbox; @@ -118,7 +121,8 @@ public: double Vex; double inAlpha; double outAlpha; - double diffusivity; + double diffusivity; + double dispersivity; }; /*! * \copydoc FvBaseLocalResidual::computeStorage @@ -295,9 +299,10 @@ public: // for thermal harmonic mean of half trans const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx); const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn); - const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn); + const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn); + const Scalar dispersivity = problem.dispersivity(globalIndexEx, globalIndexIn); - const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha, diffusivity}; + const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity}; calculateFluxes_(flux, darcy, @@ -441,7 +446,7 @@ public: static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling."); // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); - // deal with diffusion (if present) + // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity). if constexpr(enableDiffusion){ typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient; DiffusionModule::ExtensiveQuantities::update(effectiveDiffusionCoefficient, intQuantsIn, intQuantsEx); @@ -453,6 +458,19 @@ public: tmpdiffusivity, effectiveDiffusionCoefficient); + } + // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity). + if constexpr(enableDispersion){ + typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg; + DispersionModule::ExtensiveQuantities::update(normVelocityAvg, intQuantsIn, intQuantsEx); + const Scalar dispersivity = nbInfo.dispersivity; + const Scalar tmpdispersivity = dispersivity / faceArea; + DispersionModule::addDispersiveFlux(flux, + intQuantsIn.fluidState(), + intQuantsEx.fluidState(), + tmpdispersivity, + normVelocityAvg); + } // deal with micp (if present) static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling."); diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index 7d55712c7..d7f7a9ddc 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/opm/models/blackoil/blackoilmodel.hh @@ -53,6 +53,7 @@ #include #include #include "blackoildiffusionmodule.hh" +#include "blackoildispersionmodule.hh" #include #include @@ -177,6 +178,10 @@ struct EnableEnergy { static constexpr bool value template struct EnableDiffusion { static constexpr bool value = false; }; +//! disable disperison by default +template +struct EnableDispersion { static constexpr bool value = false; }; + //! by default, scale the energy equation by the inverse of the energy required to heat //! up one kg of water by 30 Kelvin. If we conserve surface volumes, this must be divided //! by the weight of one cubic meter of water. This is required to make the "dumb" linear @@ -227,7 +232,7 @@ namespace Opm { * \f] * * Since the gas and water phases are assumed to be immiscible, this - * is sufficint to calculate their density. For the formation volume + * is sufficient to calculate their density. For the formation volume * factor of the the oil phase \f$B_o\f$ determines the density of * *saturated* oil, i.e. the density of the oil phase if some gas * phase is present. @@ -288,6 +293,7 @@ private: enum { numComponents = FluidSystem::numComponents }; enum { numEq = getPropValue() }; enum { enableDiffusion = getPropValue() }; + enum { enableDispersion = getPropValue() }; static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0; static constexpr bool waterEnabled = Indices::waterEnabled; @@ -297,6 +303,7 @@ private: using PolymerModule = BlackOilPolymerModule; using EnergyModule = BlackOilEnergyModule; using DiffusionModule = BlackOilDiffusionModule; + using DispersionModule = BlackOilDispersionModule; using MICPModule = BlackOilMICPModule; public: diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh index 91fa60c31..abefea91d 100644 --- a/opm/models/common/multiphasebaseproperties.hh +++ b/opm/models/common/multiphasebaseproperties.hh @@ -80,6 +80,9 @@ struct EnableGravity { using type = UndefinedProperty; }; //! Enable diffusive fluxes? template struct EnableDiffusion { using type = UndefinedProperty; }; +//! Enable dispersive fluxes? +template +struct EnableDispersion { using type = UndefinedProperty; }; } // namespace Opm::Properties diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 54912bddf..5965442e8 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -329,6 +329,16 @@ public: return floresInfo_; } + /*! + * \brief Return constant reference to the velocityInfo. + * + * (This object is only non-empty if the DISPERC keyword is true.) + */ + const auto& getVelocityInfo() const{ + + return velocityInfo_; + } + void updateDiscretizationParameters() { updateStoredTransmissibilities(); @@ -451,6 +461,7 @@ private: Scalar outAlpha {0.}; FaceDirection dirId = FaceDirection::Unknown; Scalar diffusivity {0.}; + Scalar dispersivity {0.}; if constexpr(enableEnergy){ inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx); outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx); @@ -458,10 +469,13 @@ private: if constexpr(enableDiffusion){ diffusivity = problem_().diffusivity(myIdx, neighborIdx); } + if (simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) { + dispersivity = problem_().dispersivity(myIdx, neighborIdx); + } if (materialLawManager->hasDirectionalRelperms()) { dirId = scvf.faceDirFromDirId(); } - loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha, diffusivity}, nullptr}; + loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity}, nullptr}; } } @@ -522,15 +536,17 @@ private: jacobian_->clear(); } - // Initialize the flows and flores sparse tables + // Initialize the flows, flores, and velocity sparse tables void createFlows_() { OPM_TIMEBLOCK(createFlows); // If FLOWS/FLORES is set in any RPTRST in the schedule, then we initializate the sparse tables // For now, do the same also if any block flows are requested (TODO: only save requested cells...) + // If DISPERC is in the deck, we initialize the sparse table here as well. const bool anyFlows = simulator_().problem().eclWriter()->eclOutputModule().anyFlows(); - const bool anyFlores = simulator_().problem().eclWriter()->eclOutputModule().anyFlores(); - if ((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) { + const bool anyFlores = simulator_().problem().eclWriter()->eclOutputModule().anyFlores(); + const bool enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion(); + if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) { return; } const auto& model = model_(); @@ -539,6 +555,7 @@ private: unsigned numCells = model.numTotalDof(); std::unordered_multimap> nncIndices; std::vector loc_flinfo; + std::vector loc_vlinfo; unsigned int nncId = 0; VectorBlock flow(0.0); @@ -555,12 +572,16 @@ private: if (anyFlores) { floresInfo_.reserve(numCells, 6 * numCells); } + if (enableDispersion) { + velocityInfo_.reserve(numCells, 6 * numCells); + } for (const auto& elem : elements(gridView_())) { stencil.update(elem); for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) { unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx); loc_flinfo.resize(stencil.numDof() - 1); + loc_vlinfo.resize(stencil.numDof() - 1); for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); if (dofIdx > 0) { @@ -579,6 +600,7 @@ private: } } loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId}; + loc_vlinfo[dofIdx - 1] = VelocityInfo{flow}; } } if (anyFlows) { @@ -587,6 +609,9 @@ private: if (anyFlores) { floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end()); } + if (enableDispersion) { + velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end()); + } } } } @@ -626,7 +651,7 @@ private: // We do not call resetSystem_() here, since that will set // the full system to zero, not just our part. // Instead, that must be called before starting the linearization. - + const bool& enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion(); const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() || simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows(); const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores(); @@ -661,6 +686,11 @@ private: const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0); LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo); adres *= nbInfo.res_nbinfo.faceArea; + if (enableDispersion) { + for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { + velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea; + } + } if (enableFlows) { for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { flowsInfo_[globI][loc].flow[phaseIdx] = adres[phaseIdx].value(); @@ -823,6 +853,12 @@ private: SparseTable flowsInfo_; SparseTable floresInfo_; + struct VelocityInfo + { + VectorBlock velocity; + }; + SparseTable velocityInfo_; + using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState; struct BoundaryConditionData {