From bc03821d57f4f2ac730b2917707a7b24a03a1694 Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Wed, 25 Oct 2023 19:46:55 +0200 Subject: [PATCH] Support for mechanical dispersion --- ebos/eclalugridvanguard.hh | 4 +- ebos/eclcpgridvanguard.hh | 3 +- ebos/eclproblem.hh | 14 +++- ebos/eclproblem_properties.hh | 6 ++ ebos/ecltransmissibility.hh | 13 +++- ebos/ecltransmissibility_impl.hh | 76 ++++++++++++++++++- flow/flow_ebos_gasoil_energy.cpp | 3 + flow/flow_ebos_gasoildiffuse.cpp | 4 +- ...flow_ebos_gaswater_dissolution_diffuse.cpp | 3 + flow/flow_ebos_gaswater_energy.cpp | 5 +- opm/simulators/flow/BlackoilModelEbos.hpp | 4 + 11 files changed, 127 insertions(+), 8 deletions(-) diff --git a/ebos/eclalugridvanguard.hh b/ebos/eclalugridvanguard.hh index 31e33c270..ad699827c 100644 --- a/ebos/eclalugridvanguard.hh +++ b/ebos/eclalugridvanguard.hh @@ -189,7 +189,9 @@ public: getPropValue(), getPropValue()); + Properties::EnableDiffusion>(), + getPropValue()); // Re-ordering for ALUGrid globalTrans_->update(false, [&](unsigned int i) { return gridEquilIdxToGridIdx(i);}); } diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 0b4b89fe4..915cc144a 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/ebos/eclcpgridvanguard.hh @@ -281,7 +281,8 @@ protected: this->grid(), this->cellCentroids(), getPropValue(), - getPropValue())); + getPropValue(), + getPropValue())); globalTrans_->update(false); } diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 776f85769..a52ef1e86 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -145,6 +145,7 @@ class EclProblem : public GetPropType enum { enableTemperature = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; + enum { enableDispersion = getPropValue() }; enum { enableThermalFluxBoundaries = getPropValue() }; enum { enableApiTracking = getPropValue() }; enum { enableMICP = getPropValue() }; @@ -278,7 +279,8 @@ public: simulator.vanguard().grid(), simulator.vanguard().cellCentroids(), enableEnergy, - enableDiffusion) + enableDiffusion, + enableDispersion) , thresholdPressures_(simulator) , wellModel_(simulator) , aquiferModel_(simulator) @@ -824,6 +826,13 @@ public: return transmissibilities_.diffusivity(globalCellIn, globalCellOut); } + /*! + * give the dispersivity for a face i.e. pair. + */ + Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{ + return transmissibilities_.dispersivity(globalCellIn, globalCellOut); + } + /*! * \brief Direct access to a boundary transmissibility. */ @@ -2448,6 +2457,7 @@ private: ConditionalStorage thermalHalfTransIn; ConditionalStorage thermalHalfTransOut; ConditionalStorage diffusivity; + ConditionalStorage dispersivity; Scalar transmissibility; }; @@ -2473,6 +2483,8 @@ private: } if constexpr (enableDiffusion) *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx); + if (enableDispersion) + dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx); } }; diff --git a/ebos/eclproblem_properties.hh b/ebos/eclproblem_properties.hh index ac6cca435..4138fdd77 100644 --- a/ebos/eclproblem_properties.hh +++ b/ebos/eclproblem_properties.hh @@ -269,6 +269,12 @@ struct EnableDiffusion { static constexpr bool value = true; }; +// Enable dispersion +template +struct EnableDispersion { + static constexpr bool value = false; +}; + // only write the solutions for the report steps to disk template struct EnableWriteAllSolutions { diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index aa3219f76..76f6bfb3c 100644 --- a/ebos/ecltransmissibility.hh +++ b/ebos/ecltransmissibility.hh @@ -66,7 +66,8 @@ public: const Grid& grid, std::function(int)> centroids, bool enableEnergy, - bool enableDiffusivity); + bool enableDiffusivity, + bool enableDispersivity); /*! * \brief Return the permeability for an element. @@ -107,6 +108,11 @@ public: */ Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const; + /*! + * \brief Return the dispersivity for the intersection between two elements. + */ + Scalar dispersivity(unsigned elemIdx1, unsigned elemIdx2) const; + /*! * \brief Actually compute the transmissibility over a face as a pre-compute step. * @@ -230,6 +236,8 @@ protected: void extractPorosity_(); + void extractDispersion_(); + void computeHalfTrans_(Scalar& halfTrans, const DimVector& areaNormal, int faceIdx, // in the reference element that contains the intersection @@ -258,6 +266,7 @@ protected: std::vector permeability_; std::vector porosity_; + std::vector dispersion_; std::unordered_map trans_; const EclipseState& eclState_; const GridView& gridView_; @@ -269,8 +278,10 @@ protected: std::map, Scalar> thermalHalfTransBoundary_; bool enableEnergy_; bool enableDiffusivity_; + bool enableDispersivity_; std::unordered_map thermalHalfTrans_; //NB this is based on direction map size is ca 2*trans_ (diffusivity_) std::unordered_map diffusivity_; + std::unordered_map dispersivity_; const LookUpData lookUpData_; const LookUpCartesianData lookUpCartesianData_; diff --git a/ebos/ecltransmissibility_impl.hh b/ebos/ecltransmissibility_impl.hh index 141826bcf..317fd9659 100644 --- a/ebos/ecltransmissibility_impl.hh +++ b/ebos/ecltransmissibility_impl.hh @@ -93,7 +93,8 @@ EclTransmissibility(const EclipseState& eclState, const Grid& grid, std::function(int)> centroids, bool enableEnergy, - bool enableDiffusivity) + bool enableDiffusivity, + bool enableDispersivity) : eclState_(eclState) , gridView_(gridView) , cartMapper_(cartMapper) @@ -101,6 +102,7 @@ EclTransmissibility(const EclipseState& eclState, , centroids_(centroids) , enableEnergy_(enableEnergy) , enableDiffusivity_(enableDiffusivity) + , enableDispersivity_(enableDispersivity) , lookUpData_(gridView) , lookUpCartesianData_(gridView, cartMapper) { @@ -147,6 +149,17 @@ diffusivity(unsigned elemIdx1, unsigned elemIdx2) const } +template +Scalar EclTransmissibility:: +dispersivity(unsigned elemIdx1, unsigned elemIdx2) const +{ + if (dispersivity_.empty()) + return 0.0; + + return dispersivity_.at(isId(elemIdx1, elemIdx2)); + +} + template void EclTransmissibility:: update(bool global, const std::function& map, const bool applyNncMultregT) @@ -159,7 +172,8 @@ update(bool global, const std::function& map, const unsigned numElements = elemMapper.size(); // get the ntg values, the ntg values are modified for the cells merged with minpv const std::vector& ntg = this->lookUpData_.assignFieldPropsDoubleOnLeaf(eclState_.fieldProps(), "NTG", numElements); - const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_; + const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive(); + const bool updateDispersivity = eclState_.getSimulationConfig().rock_config().dispersion(); if (map) extractPermeability_(map); @@ -209,6 +223,13 @@ update(bool global, const std::function& map, const extractPorosity_(); } + // if dispersion is enabled, let's do the same for the "dispersivity" + if (updateDispersivity) { + dispersivity_.clear(); + dispersivity_.reserve(numElements*3*1.05); + extractDispersion_(); + } + // The MULTZ needs special case if the option is ALL // Then the smallest multiplier is applied. // Default is to apply the top and bottom multiplier @@ -314,6 +335,9 @@ update(bool global, const std::function& map, const if (updateDiffusivity) { diffusivity_[isId(elemIdx, outsideElemIdx)] = 0.0; } + if (updateDispersivity) { + dispersivity_[isId(elemIdx, outsideElemIdx)] = 0.0; + } continue; } @@ -485,6 +509,42 @@ update(bool global, const std::function& map, const diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity; } + + // update the "dispersivity half transmissibility" for the intersection + if (updateDispersivity) { + + Scalar halfDispersivity1; + Scalar halfDispersivity2; + + computeHalfDiffusivity_(halfDispersivity1, + faceAreaNormal, + distanceVector_(faceCenterInside, + intersection.indexInInside(), + elemIdx, + axisCentroids), + dispersion_[elemIdx]); + computeHalfDiffusivity_(halfDispersivity2, + faceAreaNormal, + distanceVector_(faceCenterOutside, + intersection.indexInOutside(), + outsideElemIdx, + axisCentroids), + dispersion_[outsideElemIdx]); + + applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg); + applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg); + + //TODO Add support for multipliers + Scalar dispersivity; + if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30) + // avoid division by zero + dispersivity = 0.0; + else + dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2); + + + dispersivity_[isId(elemIdx, outsideElemIdx)] = dispersivity; + } } } @@ -615,6 +675,18 @@ extractPorosity_() "(The PORO keywords are missing)"); } +template +void EclTransmissibility:: +extractDispersion_() +{ + if (!enableDispersivity_) { + throw std::runtime_error("Dispersion disabled at compile time, but the deck " + "contains the DISPERC keyword."); + } + const auto& fp = eclState_.fieldProps(); + dispersion_ = fp.get_double("DISPERC"); +} + template void EclTransmissibility:: removeSmallNonCartesianTransmissibilities_() diff --git a/flow/flow_ebos_gasoil_energy.cpp b/flow/flow_ebos_gasoil_energy.cpp index 5e70d9a77..74af830db 100644 --- a/flow/flow_ebos_gasoil_energy.cpp +++ b/flow/flow_ebos_gasoil_energy.cpp @@ -70,6 +70,9 @@ struct LocalResidual { using type = B template struct EnableDiffusion { static constexpr bool value = true; }; +template +struct EnableDispersion { static constexpr bool value = true; }; + }} namespace Opm { diff --git a/flow/flow_ebos_gasoildiffuse.cpp b/flow/flow_ebos_gasoildiffuse.cpp index 31b2fe070..ae62a01c2 100644 --- a/flow/flow_ebos_gasoildiffuse.cpp +++ b/flow/flow_ebos_gasoildiffuse.cpp @@ -41,10 +41,12 @@ struct Linearizer { using type = Tpf template struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; - template struct EnableDiffusion { static constexpr bool value = true; }; +template +struct EnableDispersion { static constexpr bool value = true; }; + //! The indices required by the model template struct Indices diff --git a/flow/flow_ebos_gaswater_dissolution_diffuse.cpp b/flow/flow_ebos_gaswater_dissolution_diffuse.cpp index 16626b094..64b651277 100644 --- a/flow/flow_ebos_gaswater_dissolution_diffuse.cpp +++ b/flow/flow_ebos_gaswater_dissolution_diffuse.cpp @@ -48,6 +48,9 @@ struct LocalResidual { template struct EnableDiffusion { static constexpr bool value = true; }; +template +struct EnableDispersion { static constexpr bool value = true; }; + template struct EnableDisgasInWater { static constexpr bool value = true; diff --git a/flow/flow_ebos_gaswater_energy.cpp b/flow/flow_ebos_gaswater_energy.cpp index d7028682f..da363133f 100644 --- a/flow/flow_ebos_gaswater_energy.cpp +++ b/flow/flow_ebos_gaswater_energy.cpp @@ -46,7 +46,10 @@ template struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; template -struct EnableDiffusion { static constexpr bool value = false; }; +struct EnableDiffusion { static constexpr bool value = true; }; + +template +struct EnableDispersion { static constexpr bool value = true; }; template struct EnableEnergy { diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 45d73d699..70dc5147d 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -133,6 +133,10 @@ template struct EnableMICP { static constexpr bool value = false; }; +template +struct EnableDispersion { + static constexpr bool value = false; +}; template struct EclWellModel {